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(]j Abstract 

The subject of this study is the long-time equilibration dynamics of a strongly disor- 

^ dered one-dimensional chain of coupled weakly anharmonic classical oscillators. It is 

shown that chaos in this system has a very particular spatial structure: it can be viewed 

r^ as a dilute gas of chaotic spots. Each chaotic spot corresponds to a stochastic pump 

Ch which drives the Arnold diffusion of the oscillators surrounding it, thus leading to their 

tXi relaxation and thermalization. The most important mechanism of relaxation at long dis- 

^ tances is provided by random migration of the chaotic spots along the chain, which bears 

analogy with variable-range hopping of electrons in strongly disordered solids. The cor- 

resonding macroscopic transport equations are obtained. 
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cn 1. Introduction 
> 

^^ Anderson localization [T| is a general phenomenon occurring in many linear wave- 

Q like systems subject to a disordered background (see Ref. [2| for a recent review). It 

l/~) is especially pronounced in one-dimensional systems, where even an arbitrarily weak 

I/-) disorder localizes all normal modes of the system ll3]|4|. Anderson localization in linear 

O systems (single-particle problems in quantum mechanics) has been thoroughly studied 

^ over the last 50 years; its physical picture is quite clear by now (although some open 

r| questions still remain), and even rigorous mathematical results have been estabhshed. 

• ^H The situation is much less clear in the presence of nonlinearities/interactions. 
rN One of the simplest systems where the effect of a classical nonlinearity on the An- 

C^ derson localization can be studied, is the disordered nonlinear Schrodinger equation 
(DNLSE) in one dimension, discrete or continuous. It has attracted a lot of interest in 
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the last few years, including both stationary and non-stationary problems, as well as the 
problem of the dynamic stability of stationary solutions I5l|6l|7l[8l|9l[10l[ni[l2l[l3l[ll 

IBl [IS 113 [II [13 lia ini 122 El El 123 |2i 123 IM 113 111 E^ 

This interest was especially motivated by the experimental observation of the Ander- 
son localization of light in disordered photonic lattices |[38l [39l |40|| . where DNLSE 
describes light propagation in the paraxial approximation, and of cold atoms in disor- 
dered optical lattices Pni42ll43l . for which DNLSE provides the mean-field descrip- 
tion. Other systems have also been studied, such as a chain of anharmonic oscilla- 
tors ||3|I5llMlll5ll20ll2ll[32l[3l[35l|46l, a chain of classical spins |[51|43, a nonlinear 
Stark ladder mijia. 

In contrast to linear problems which can always be reduced to finding the normal 
modes and the corresponding eigenvalues, the problem of localization in a nonlinear 
system can be stated in different inequivalent ways. For example, one can study solutions 
of the stationary nonlinear Schrodinger equation with disorder llT2l [T4l [28l ; however, in 
a nonlinear system they are not directly related to the dynamics. Another possible setting 
is the system subject to an external perturbation or probe; studies of the transmission of 
a finite-length sample where an external flow is imposed |6, 8) |3 [IHl [IH [13 123 l27l . 
or dipolar oscillations in an external trap potential [18|, fall into this category. Most 
attention has been paid to the problem of spreading of an initially localized wave packet 
of a finite norm |[3[ai3[l3l[l5l[l6l[ni|50l|20l|2l]|221|23l|43|43l2a|23|^[3ll 
l32l l33l |34]| (in hnear systems it remains exponentially suppressed at long distances for 
all times). This setting corresponds directly to experiments 1(381 [39l l40l |4T1 l42l . A 
closely related but not equivalent problem is that of thermalization of an initially out- 
of-equilibrium system (in a linear localized system thermalization does not occur). The 
difference between these two settings is that in the latter the total energy stored in the 
system is proportional to its (infinite) size, while in the former the infinite system initially 
receives a finite amount of energy. The problem of thermalization has also received 
attention Il44l l45l l47l |26| . and it is the main subject of the present work (although the 
problem of wave packet spreading will also be briefly discussed). Still, despite a large 
body of work, detailed understanding of the equilibration dynamics in one-dimensional 
disordered nonlinear systems is still lacking. 

On the one hand, the direct numerical integration of the differential equation shows 
that an initially localized wave packet of a finite norm does spread indefinitely, its size 
growing with time as a sub-diffusive power law |13 13 [13 [I6ll20l I2I1|321 ISS |33|. Ex- 
istence of different regimes of spreading, depending on whether the system is in the 
regime of strong or weak chaos, has been discussed |[33ll35]| . Several authors suggested 
a macroscopic description of the long-time dynamics in terms of a nonlinear diffusion- 
type equation resulting in sub-diffusive spreading of the wave packet fT^ I^Tl [^ ; how- 
ever, arguments used to justify this equations are based on an a priori assumption of 
spatially uniform chaos. 

On the other hand, it was argued that among different initial conditions with a finite 



norm, at least some should exhibit regular quasi-periodic dynamics, and thus would not 
spread indefinitely |l5l[T3l|37j|. For finite-size systems, scaling of the probability for the 
system to be in the chaotic or regular regime has been studied numerically |36|. When 
scaling arguments are applied to the results of the numerical integration of the equation 
of motion, they indicate slowing down of the power-law spreading |[34l . Rigorous math- 
ematical arguments support the conclusion that at long times the wave packet should 
spread (if spread at all) slower than any power law [50| . 

The aim of the present work is the analysis of statistical properties of weak chaos in 
the discrete one-dimensional disordered nonlinear Schrodinger equation with the initial 
conditions corresponding to finite norm and energy densities. The main focus is on long- 
time dynamics and relaxation at large distances, when no stationary superfluid flow can 
exist in the system lITTl l27l and the dynamics is chaotic. It is often assumed that upon 
thermalization chaos has no spatial structure, and all sites of the chain are more or less 
equally chaotic; here it is argued not to be the case{^ Namely, in the regime of strong 
disorder and weak nonlinearity chaos is concentrated on a small number of rare chaotic 
spots (essentially the same picture was also proposed in Ref. Il47l for a disordered chain 
of coupled classical spins, but it was not put on a quantitative basis). A chaotic spot 
is a collection of resonantly coupled oscillators, in which one can separate a collective 
degree of freedom (namely, their relative phase whose dynamics is slow because of the 
resonance), which performs chaotic motion. Under the conditions of weak coupling be- 
tween neighboring oscillators and weak nonlinearity (i. e., Anderson localization being 
the strongest effect), the distance between different chaotic spots is much greater than 
the typical size of the spot. The stochastic motion of the collective degree of freedom of 
each chaotic spot acts as a stochastic pump, i. e. drives the exchange of energy between 
other oscillators, non-resonantly coupled to the spot, which corresponds to the Arnold 
diffusion. This represents the main mechanism for relaxation and thermalization of the 
oscillators, as well as for the transport of conserved quantities (energy and norm). An 
important role is played by the fact that chaotic spots can migrate along the chain, as the 
collective degree of freedom may get in and out of the chaotic region of its phase space. 
A similar phenomenon has recently been proven to exist in a chain of weakly coupled 
pendula [51 1. This migration also bears analogy with the variable -range hopping of elec- 
trons in strongly disordered solids |[52ll53l . Still, there is an important difference that the 
chaotic spot does not carry any energy or norm; it only drives relaxation of oscillators 
surrounding it. 

The main technical challenge in this work is the analysis of high orders of the per- 
turbation theory, which is necessary both to separate the collective degree of freedom 
performing the chaotic motion, and to couple other oscillators to this degree of freedom. 
The perturbation theory is divergent, as chaos is a non-perturbative phenomenon. The 
divergence occurs because of overlapping resonances, according to Chirikov's criterion 



' We will call the motion of an oscillator the more chaotic the faster it loses the phase memory. 
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of chaos. In this work, the probabiUty for resonances to occur is estimated in each order 
of perturbation theory. The subsequent treatment of each resonance and description of 
the associated chaotic motion is based on the stochastic pump model of the Arnold dif- 
fusion, which was proposed for systems with few degrees of freedom |l54l|55l. Thus, the 
present work is not more rigorous than the stochastic pump model; it does not add any- 
thing new to the current understanding of how chaos is generated in nonlinear systems; 
it rather deals with the statistics of chaos in a spatially extended system with an infinite 
number of degrees of freedom and local coupling between them. 

The main quantitative result of the present work is the system of macroscopic equa- 
tions which describe transport of the conserved quantities (energy and norm) along the 
chain. Such tansport determines equilibration of the chain at large distances, which oc- 
curs at long time scales. These equations are of nonlinear-diffusion type, and explicit 
expressions for the transport coefficients are derived. The macroscopic equations are 
valid at distances exceeding a certain length scale, which is also explicitly estimated. 

The paper is organized as follows. In Sec. |2] the model is formulated and the main 
assumptions are discussed. Sec. [3] summarizes the main results. In Sec. |4] we quali- 
tatively describe the physical picture and the main ingredients of the solution; various 
implications and some related issues are discussed. In Sec. [5] we analyze the chaotic 
phase space in few-oscillator configurations. In Sec.[6]perturbation theory is developed 
and statistics of different terms is analyzed. In Sec.[7]this perturbation theory is used to 
study the statistics of the chaotic phase volume for an arbitrary number of oscillators. 
Sec. [8] is dedicated to the analysis of the coupling between chaotic spots and other os- 
cillators, and to the description of Arnold diffusion. In Sec. |9]the macroscopic transport 
coefficients are found. 

2. Model and assumptions 

The main subject of this study is the one-dimensional discrete nonlinear Schrodinger 
equation with diagonal disorder: 

dlj/n H, 2 

i—7- =<J}n^n-Tl^{^n-\ + >JJn+\) + g^n^n- (2-1) 

Here the integer n runs from -oo to +00 and labels sites of a one-dimensional lattice. To 
each site a pair of complex conjugate variables i/^,,, ip*^ is associated. 



Eq. ( 2. 1 1 together with its complex conjugate are the Hamilton equations, 

dipn _ OH d{i4tl) _ dH 
dt d{i>Jj*„y dt dipn 

corresponding to the classical Hamiltonian 

//({/(A*}, {lAn}) = ^ loJnK^n - tA (lA^^+l + K+l^n) + X ('f'nf^l 



(2.2) 



(2.3) 



where iij/*^ is the canonical momentum conjugate to the coordinate iJ/„. For this i/^„ must 
have a dimensionality of (action) ^^^. One can pass to another set of canonical variables, 
the actions /„ and the phases 0,,, as 

<An = V^^"''^". (2-4) 

so that the Hamiltonian becomes 

H{{In}, {(f>n}) = ^ lojJr, - 2tA V44+1 COS(0„ - (^„+i) + ^^l ^ 



(2.5) 



and the equations of motion are 



d(pn _ dH dl„ _ dH 
dt din ' dt d(pn 



(2.6) 



It is convenient to introduce the action-dependent frequency of each oscillator, 

ibniln) ^(^n+ gh, (2.7) 

to be distinguished from the bare frequencies 6t»„. 



The first term on the right-hand side of Eq. (2.1 1 represents the diagonal disorder. 
The bare frequencies w„ are assumed to be independent random variables uniformly 

distributed over the interval 

A A 



The second term on the right-hand side of Eq. (2.1 ) will be referred to as tunnelingjj It 



is convenient to measure its strength with respect to the disorder, so the dimensionless 
parameter t is introduced. Together with the disorder term, the tunneling term constitutes 
the linear part of the problem, whose normal modes are determined from the eigenvalue 
equation 

(a» - W„)(A„ + T^{ipn-\ + lAn+l) = 0. (2.9) 

As we are dealing with a one-dimensional system, any eigenfunction has an envelope, 
exponentially localized around some site h, if/^ ~ g-^l^-wl/^ x^g smaller r, the stronger 
the localization. We assume 

T«l. (2.10) 

In this case k x ln(l/2eT) » ijjso each eigenstate is essentially localized on a single 
site, with weak perturbative tails on the neighboring sites. 



"May the reader not be confused by this "quantum-mechanical" term: here we are always dealing with a 
classical Hamiltonian system. This term is the standard one in studies of Anderson localization on a lattice. 

^^This estimate can be obtained, e. g. from the Thouless relation between the localization length and the 
density of states |56| , approximating the latter by the box of the width A. 
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The third term on the right-hand side of Eq. ( |2.1[ ) represents the nonlinearity whose 
strength is governed by the coupling constant g > OJj As mentioned above, we are 
interested in the case when the nonlinearity is weak. However, it would be meaningless 



to simply say that g is small. Indeed, Eq. ( 2. 1 1 is invariant under the change 



tpn^CljJn, S ^ -^, (2-11) 

where C is an arbitrary constant. Thus, it is the product g\tj/n\^, the nonlinear frequency 
shift of the nth oscillator, which can be compared to other frequency scales. This non- 
linear shift is determined by g, as well as by the initial conditions for Eq. ( |2. 1[ ). We will 
assume that it is much smaller than the disorder and denote the corresponding dimen- 
sionless small parameter by p: 

^^~p«l. (2.12) 



Eq. (2. 12 1 also expresses another assumption of the present work, that the typical 
scale of li/'np is about the same for all oscillators; in other words, there is no significant 
correlation between \ifr„\^ and a>„. As an example of the opposite one could prepare an 
initial condition where most of the norm \ifr„\^ is put on the oscillators with a»„'s near 
the bottom of the band. This case is more complicated (e. g., it includes all physics 
related to superfluidity, such as Bose condensate being the ground state of the classical 



Hamiltonian (2.3 1, as well as large localization length for low-frequency phonons |57| ); 



it is beyond the scope of the present work. The criterion can be conveniently formulated 



in terms of the two integrals of motion of Eq. (2. 1 1: the total energy, given by Eq. (2.3 1, 
and the total action (norm). 



£ m 



(2.13) 



In the situation favoring oscillators with frequencies near the bottom of the band one 
obtains H ~ -{A/2)Itot, since the dominant contribution to both H and It„t comes from 
the oscillators with a>n close to -A/2 (the coupling and the nonlinear terms in H are 
assumed to be small). In the case with li/^np independent of co„ there is a significant 



cancellation between contributions a)n\tfr„\^ to the Hamiltonian (2.3 1 with cj,, > and 
cOn < 0, SO assuming all oscillators to be excited on the same scale is equivalent to 
assuming \H\/Itot ^ A/2. In the following, we will assume the strong inequality, 

^ « A, (2.14) 

'tot 



''Generally speaking, the sign of g is not important for the dynamics, since the change g — > -g is 
equivalent to i^„ — » {-l)"if/l, a>„ — » -a»„; the latter has no effect since the disorder distribution is symmetric. 
However, later we will be interested in the thermodynamics of the chain, so it is convenient to have the 
Hamiltonian bounded from below. 
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which is made for pure technical convenience only. It does not affect the qualitative 
conclusions of the work, but simplifies the calculations since it effectively decouples 
different realizations of disorder and different initial conditions. 
Finally, we will assume 

H<0. (2.15) 

It turns out that for certain values of H and Itot the system cannot thermalize, which can 
be seen from very basic thermodynamic arguments (see Ref. [58 1 and [Appendix A for 



details). Condition ( 2. 15 1 in combination with Eqs. (2. 12 1, ( 2. 14 ) ensures that the system 



is not in that non-thermal region. 

When all three terms on the right-hand side of Eq. (2.11 are present, the equation 
has only two integrals of motion: the total energy, H, and the total action, Itot - Z„ ^«- 
However, if r = or g = 0, the system is integrable. Thus, one is free to choose any 
of the two as the perturbation destroying the integrability. The most natural choice is 



the tunneling. Then, each /„ being an integral of motion at t = 0, Hamiltonian (2.5 1 
is already written in the correct variables. It may seem that this implies the condition 
T <^ p (since p is the effective strength of the nonlinearity). Indeed, as will be seen below 
(Sec. |5]l, some elements of the physical picture admit a simpler quantitative description 
at T «; p. However, the main results of this paper are obtained using the perturbation 
theory both in t and g, and are valid for any relation between t and p, as long as both 



are small (on the issue of smallness, see also Sec. 4.5.6) 



3. Main results 

The first statement is that for general initial conditions, satisfying the assumptions of 
the previous section, the system locally thermalizes with a finite relaxation time. Local 
thermalization means that on a finite, sufficiently long segment of the chain [the corre- 
sponding length scale L* is given explicitly below, Eq. ( |3.7| )], the actions and the phases 
are distributed according to the grand canonical distribution, e~^'^^''^'""\ To respect the 
two conserved quantities, the total energy H and action Itot of the chosen segment, two 
Lagrange multipliers are introduced: yS = 1/r is the inverse temperature, and // has the 
meaning of the chemical potential, although it has the dimensionality of a frequency. 
The issue of thermalization is discussed in more detail in Sec. 14.5.21 

The values of ju and T, which define the thermal distribution on each segment, are 
determined by the total norm and energy initially contained on this segment, that is, by 
the initial conditions. It turns out that for some values of energy and norm (namely, 
when the energy is too high for a given norm), the corresponding ju and T do not exist. 



as discussed in Ref. [,58il and Appendix A so the system may equilibrate in the micro- 



canonical sense, but still not be thermal. The assumption of negative total energy, made 



in the previous section [Eq. (2.151], ensures that such problem does not arise, and the 
system indeed thermalizes. However, even when the local thermalization has occurred, 
generally speaking, // and T still have different values for different segments, depending 
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on the initial conditions, so the system is still not in the global equilibrium. The results 
described in the following, concern precisely the global equilibration of the system. 

To account for the dependence of ju and T on the position, we introduce a continuous 
variable, n ^ x, and assume a dependence on the position, T = T{x), /u = fiix), smooth 
on the scale L* . Let us define the smoothed macroscopic energy and action densities. 



X+LJ2 



* n=x-L.I2 



n=x-L,/2 



(3.1a) 
(3.1b) 



At each point x they are related to the local values of T{x),fj.{x) by the standard thermo- 



dynamics (see Appendix A i. It is convenient to define dimensionless variables 



P = 



gl 



h = 



gn± 

A2' 



(3.2) 



Under the assumptions of the previous section, namely, T«;l,p«;l, |/j|«:p, the 
thermodynamic relations take a very simple form (the condition \h\ -« p turns out to be 
equivalent to // < 0, |/i| » A): 



l)U|A' 



h 









(3.3) 



In the global equilibrium, T and ju must be constant over the entire chain. As the 
values ot/uix) and T{x) are determined by the densities I{x) and 'K(x), during the global 
equilibration norm and energy must be transported along the chain. Thus, the corre- 
sponding currents, J-^ and f^, must be flowing during the equilibration. The conserva- 
tion laws for the total energy and norm impose the corresponding continuity equations: 



d_ 
dt 



I 



d_ 
dx 



J' 



(3.4) 



The currents must vanish for T{x),iJ.{x) = const. Thus, they can be expanded in the 
gradients, and only the first term is retained: 



J^ 
J^ 



-(T{T,ii) 



d 
dx 



-1/r 



(3.5) 



cr(T,^) is a 2 X 2 matrix of the transport coefficients which will be called conductiv- 
ity. Since -jijT and 1/r are thermodynamically conjugate to I and "K, the matrix cr is 
symmetric by virtue of the Onsager relations. The statement about the finiteness of cr, 
equivalent to the validity of the gradient expansion, and to the normal character of diffu- 
sion in the system, is already not trivial. Indeed, cr could be zero or infinite, leading to 
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anomalous diffusion; in Ref. fTP\ cr = was suggested. Our result is that cr is finite, and 
given by 



A3 

Cr(r,yU) 



1 0(pA) 
0(pA) 0(A2) 



exp(-C'(!=<iM],„2_Li„i\. (3^6a) 



It depends on T and fi only through the combination p = ^r/(|ju|A), which is the con- 



sequence of the condition \iu\ » A, or, equivalently, of Eq. (2.141. The argument of the 
exponential is written with the logarithmic accuracy. The function C'(ln(l/p)/ ln(l/T)) 
is bounded from above and below. 



1/3 



[l+ln(l+x)]-^ 



^ C'(x) ^ 



[l+ln(l+x)] 



2" 



(3.6b) 



Finally, the power p is bounded by 



^ 7? ^ 3. 



(3.6c) 



The author has not been able to calculate C'{x) and p explicitly. Given this uncertainty, 
it would make little sense to calculate slow prefactors. The prefactor A^/g^ is put to 
keep the correct dimensionality. Also, an order-of-magnitude estimate for the relative 
magnitude of different matrix elements of cr could be established. 



Equations (3.2 1 -(3.6c i form a closed system of macroscopic equations describing 
action and energy transport along the chain and constitute the main quantitative result of 



the present work. Eq. ( 3.6a i is already sufficient to conclude that the dependence of cr on 
the integrability-breaking parameters r and p is stronger than any power law, but weaker 
than a stretched exponential which would be the expected dependence of the Arnold 
diffusion rate on the integrability-breaking perturbation Il54ll55ll59l . This phenomenon 
is known as fast Arnold diffusion [|60il . and arises because of a large number of degrees 
of freedom in the system. 



The above expression for cr{T,p), Eq. (3.6a I, does not imply averaging over many 
disorder realizations. In fact, its inverse, cr^\ is a self-averaging quantity. The length 
scale Lt at which this self-averaging occurs can be viewed as boundary separating macro- 
scopic and microscopic distances. It can be estimated as 



InL, -C'ln^ 



1 

TPp' 



(3.7) 



Segments much longer than L* contain a sufficient number of sites to effectively replace 
the average over realizations. Transport properties of a segment shorter than L* are 
strongly realization-dependent, and fluctuations are of the same order as the average. 
This situation is analogous to that with hopping conduction in one dimension [61]. 

The macroscopic transport equations given above describe equilibration of a chain 
whose average action I and energy 'H per oscillator are finite. Still, having written 
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these macroscopic equations, one can always prepare an initial condition J(x), T^Cx) in 
the form of a finite-size packet, and see how it spreads according to the equations. In 



Appendix B the equations are shown to give the size of the packet, L{t), to grow as 



L(,).^[(C'r'lnfl'/3_ (3_g^ 



This is slower than any power law, in agreement with Ref. IdOII . and in disagreement 
with Refs. |I71|9l[15l[^|20l|2T][32l[3l[35l. Possible reasons for this disagreement 



will be discussed in Sec. 4.5.6 Note, however, that as the packet size L{t) increases, the 



corresponding density inside the packet, p oc l/L{t), decreases. Thus, the length scale L*, 



given by Eq. (3.7 1, very quickly diverges so after some time we will inevitably have 
L{t) < Lh.. Since the macroscopic equations are only valid at distances exceeding L*, 
they lose their applicability as the packet spreads. 

4. Qualitative picture 

4. 1. Search for chaos 

Appearance of chaos in systems with a few degrees of freedom has been thoroughly 
studied in the past [54, 55]|59l; the main steps will be reviewed in Sec. [5] The essential 
ingredients are the following. 

1 . When an integer linear combination of several frequencies a)„(/„) of the integrable 
system vanishes for some values of 7„, that is mia),^ + . . . + niMOJu^ = (guiding 
resonance, in Chirikov's terminology f54l), one can separate a slow degree of 
freedom 0, which is the corresponding linear combination of the phases: = 
m 1 0„ , + . . . + m/v0„^ . Exactly at resonance this degree of freedom has no dynamics 
at all (in the unperturbed system). 

2. If the perturbation has a term resonant with 0, the latter acquires a non-trivial 
dynamics whose characteristic frequency Q. is proportional to the square root of 
the perturbation strength. The corresponding phase space has a separatrix. 

3. Chaos arises upon destruction of this separatrix by another term of the pertur- 
bation which oscillates at some characteristic frequency ru. The volume of the 
resulting chaotic region around the separatrix in the phase space of the slow de- 
gree of freedom (the so-called stochastic layer) (i) is proportional to the strength of 
the separatrix-destroying perturbation (up to a logarithmic factor), and (ii) is expo- 
nentially small, ~ e'^^^, when tir » O (which is typically the case). Destruction 



of the separatrix will be also reviewed in more detail in Sec. 5.2 



In the particular case of Eq. (2.1 1 the perturbation conserves the total action and energy, 
so the minimal number of oscillators, necessary to generate chaos, is three. 

To determine the main contribution to the chaotic phase volume, we need (i) to iden- 
tify the guiding resonances which have the largest resonant perturbation strength and 
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thus the highest frequency Q.; (ii) to find the separatrix-destroying terms of the pertur- 
bation which oscillate at the smallest possible frequency vr and are not too weak (layer 
resonances, in Chirikov's terminology ll54l ). It is important that the mentioned pertur- 



bation terms include not only those directly appearing in the Hamiltonian (2.3 1 or (2.5 1, 
but also terms of the kind 

V^'K (lA*^ . . . (A^^iA/i, ■ ■ ■ (Aha. + "As, • ■ • 'Af,„'A«i • - - 'A«;v) 

= 2% yjlny ...InN^i ■■■Ihn COS (0,,, -I- . . . -I- 0„^ - (^/j, - ... - (f>fi^) . 

Such ' W-particle" terms are effectively generated in high orders of the perturbation the- 
ory in T and g, which is developed in Sec. [6] They oscillate at the frequency m = 
cbni -I- ... -I- a)„^ - a)fn - ... - cbfiN^ which is easier to make small for larger A^, as dis- 



cussed in Sec. 6.6 However, the corresponding strength 7C contains a high power of the 
coupling constants. 

Consider the simplest resonance, that of two oscillators: 

iOn + gin ^ O^n' + gin' ■ (4.2) 

Since the typical values of actions /„, /„/ ~ r/|yu|, this equation can be satisfied only if 

gT 

\a)n-a)n'\ ~ 7T- (^.3) 



Two arbitrary sites have \a>„-CL),f | ~ A, so for weak nonlinearity, gT/\iu\ «; A [Eq. ( 2. 12 1], 
one has to look for special pairs. If we fix n, the appropriate partner n' will be found 
typically at distances \n - n'\ ~ \^\KI{gT) = I /p. This distance gives also the minimal 
order of the perturbation theory in the tunneling which couples two such oscillators 
and produces a separatrix, so that O x Vr^. The thickness of the corresponding 
stochastic layer, oc e~^^^, will be very small. The same occurs if one tries an A/^-oscillator 
resonance: even though the equation 

miiconi +gl„i) + ■■■+ mNioj„„ + gln^) = (4.4) 



is easier to satisfy than Eq. (4.2 1, a high order of the perturbation theory is required to 
couple these oscillators altogether, and thus Q. is again proportional to a high power of t, 
which then enters as e~^/^. 

The "cheapest" way to produce a moderately small Q. is to find a pair of nearest 
neighbors, n' = n ± I, such that \a>n - cOn+\\ ~ gT/\fJ-\- Such pairs are rare, they occur 
every 1 /p sites on the average, but they have Q oc ^r. The low density of these pairs 
gives a small factor ~ p in the total chaotic phase volume, but it is much less dramatic 
than T^^f' exponentiated. It is here that the assumption t «; p is necessary. Indeed, in 
the opposite case one would have to solve the linear part of the problem first and find 
the two eigenfrequencies, whose difference cannot be smaller than tA (the so-called 
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level repulsion). These eigenfrequencies should then be substituted into Eq. (4.3 1 which 
cannot be satisfied if t » p. In this latter case the resonant pair should be made out 
of oscillators separated by some distance /, so that their effective coupling is ~ t' < p, 



which can always be satisfied for a sufficiently large / (see also Sec. 6.7 1. 

Next, one has to identify the separatrix-destroying perturbation which produces the 
thickest stochastic layer. If we simply consider coupling between one of the sites of 
the resonant pair n,n + I and another oscillator n', the frequency of this perturbation is 
m = \ajn - ci>n'\- Typically, this frequency is ~ A, so the effect of the perturbation is 
exponentially small. Still, out of many available oscillators, it is possible to find one 
with a sufficiently close frequency, so that nr is also small. However, this oscillator 
will be far away, so the stochastic layer thickness will be proportional to a high power 
of T. In high orders of perturbation theory many-oscillators terms can also be generated. 
Thus, one has to find an optimum between decreasing m in e~^^^ and multiplying it 
by a large power of r. An analogous optimization procedure has been considered by 
Chirikov and Vecheslavov [ 601 . For our problem the resulting width of the stochastic 
layer is proportional to e"'" ^^/'^\ 

Again, the way to bypass this smallness is to notice that the above arguments concern 
a resonant pair in a typical surrounding. There is, however, a small fraction of pairs 
which have a neighbor with a close frequency, i. e. resonant triples. Clearly, it makes 
no sense to tune the frequency of the third oscillator more precisely than the frequency 
mismatch of the resonant pair, so the density of such resonant triples is ~ p^. Moreover, 
when all three frequencies are close, the chaotic region of the phase space does not have 
to be just a thin layer around the destroyed separatrix, but can be of the size of the 
separatrix itself. Thus, the chaotic region of the {h,l2,h} space is concentrated around 
the straight line a»i +gh - (^2+gh = (^3 +gh- A resonant triple will be called a chaotic 
spot if its dynamical variables fall in the chaotic region. The probability of this latter 
event, which turns out to be ~ r/p if t «; p and ~ 1 if r » p, multiplied by the density 
of resonant triples, ~ p^, gives the density of the chaotic spots, w ~ min{Tp,p^}. This is 
discussed in detail for t «: p in Sec.|5.4[ 



4.2. Arnold diffusion 

The key property of a chaotic spot is that the Fourier transform of i/^„(?) has a con- 



tinuous frequency spectrum, as discussed in Sec. 8.1 (see also Sec. 4.5.11. When the 
stochastic layer around the separatrix is thin, the spectrum can be evaluated explicitly 
via Melnikov- Arnold integral 1 62 |63l |54 1 . The spectrum is peaked around a)„ with 
the characteristic width Q., and has an exponential tail, g-l'^-'^«l/^ away from a)„. This 
should be contrasted with the case of an integrable system where the Fourier spectrum 
of the motion is made of 5-peaks with zero width. 

The chaotic spot exerts a random force on other oscillators, thereby allowing them 
to exchange action and energy with the spot and among themselves. This exchange 
occurs by small random increments, so the relevant picture is that of random diffusion 
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in the multi-dimensional space of actions {/„). This corresponds to the standard picture 
of stochastic pump driving Arnold diffusion |[54l l55l . This diffusion (i) must respect 
the total action and energy conservation laws, and (ii) cannot "turn off" the chaotic 
spot, since the dynamic variables of the spot oscillators cannot cross the stochastic layer 
boundary |54|, so that the system stays on the guiding resonance 2^„m„6J„ = 0. The 



corresponding diffusion equation is analyzed in Sees. 8.2 and 8.3 

To study the diffusion of an oscillator, located at large distance from the nearest 
chaotic spot, one should find a suitable perturbation term of the form ( |4.1| ) involving 
both this oscillator and an oscillator belonging to the spot (and, typically, many others). 
Again, one seeks (i) to be close to resonance with the spot oscillation frequency in order 
to avoid suppression e~^^^, and (ii) to avoid too high orders of the perturbation theory, 
since the diffusion coefficient is proportional to the perturbation strength squared. Still, 
if the distance between the chaotic spot and the oscillator is L, the smallness r^^ in the 



diffusion coefficient is already guaranteed. This is discussed in Sec. 8.4 
For a given number N of oscillators, the conditions 

N N N 

Y,h=NI, 2](w„/„ + |/2j = ^'K, 2 m„(a;„ + g/„) = 0, /i, . . . ,/iv > 0, 

n=l n=\ n=l 

(4.5) 
define a region on an (A'^ - 3) -dimensional hypersurface whose {N - 3)-dimensional 
volume is finite. Hence, it takes a finite time for the system to explore this region and 
to establish the micro-canonical distribution. This time, however, grows exponentially 
with N, as long as only one chaotic spot is present. 

For a longer segment, such that it hosts two chaotic spots, the oscillators which 
are close to one spot equilibrate among themselves, and those close to the other spot - 
among themselves, much faster than equilibration between the two groups of oscillators 
occurs. Indeed, in order to transfer action and energy between the two groups, one has to 
find a perturbation term which couples oscillators from both groups to one of the spots. 
Since the coupling decays exponentially with distance, the best way to do this would 
be to couple oscillators which are near the boundary between the two groups, i. e., in 
the middle between the two spots. As mentioned in the previous subsection, the typical 
distance between chaotic spots is ~ 1/ min{ Tp,p^), so the characteristic time scale of 
equilibration is of the order of (i/j)i/™nf^P.P"l^ which is enormously long. For the long 
chain with many spots it may become even worse if one considers anomalously long 
segments of the chain, not containing chaotic spots: even though such segments ai^e rare, 
they efficiently block the transport along the chain. 

The above picture assumes the positions of chaotic spots to be fixed. There is, how- 
ever, a mechanism which provides a faster way of equilibration at long distances. Con- 
sider first a chaotic spot and the closest resonant triple which is not chaotic, assuming 



T <^ p. As discussed in Sec. 4.1 the typical distance to this triple, ~ 1/p^, is smaller 
than the distance to the next chaotic spot, ~ l/(Tp). Thus, the diffusion coefficient for 
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the actions of the triple is small as t^^p . As the triple explores its phase space in the 
course of Arnold diffusion, it may hit its own chaotic region. Importantly, this requires 
the same amount of time as the equilibration itself jj Thus, a new chaotic spot is created, 
relatively close to the old one. However, now each of them can diffuse out of its stochas- 
tic layer, driven by the other one. If it happens to the old one, this means that the system 
has effectively switched from one guiding resonance to another one at the point of their 
intersection in the phase space. This switching is also accompanied by a jump of the 
chaotic spot in the real space. 

It is important to realize that this jump does not correspond to a transfer of any action 
or energy: the chaotic spot does not carry any action or energy by itself; these are stored 
in the numerous oscillators surrounding the spot and transferred from one oscillator to 
another by Arnold diffusion, while the chaotic spot only drives this diffusion. Thus, 
to transfer action and energy between two groups of oscillators, which already have 
equilibrated among themselves, it is quicker to wait until the chaotic spot moves closer 
to the boundary between the two groups and drives the Arnold diffusion there, rather 
than doing it by diffusion driven from a remote spot. 

If now the consideration is extended on triples which are not well in resonance, but 
have a noticeable frequency mismatch, this mismatch will translate into an activation 
barrier E, so the time required for a triple to reach its chaotic region will be longer than 
the equilibration time by a factor e^^^ . However, since the typical activation barrier of 
an arbitrary triple is ~ \fj.\A/g, hopping to a typical nearest neighbor will only take the 
time e^/P, much shorter than (1/t)^^^ . This means that resonant triples, on the one hand, 
host chaotic spots for most of the time, but on the other, are not the ones that determine 



the transport. Thus, one has to return to the beginning of Sec. 4.1 and look for guiding 
and layer resonances in a typical environment. For this, many-oscillator resonances have 
to be considered, and thus, high orders of the perturbation theory have to be analyzed. 

4.3. Variable-range hopping of chaotic spots 

Even when the guiding resonance involves several oscillators, it is still convenient 
to assign it to a single site of the chain (defined, e. g., as the leftmost site participating 
in the resonance). Then, each site n can be characterized by w„, the chaotic fraction 
of the thermally-weighted phase volume, summed over all possible guiding resonances 



involving the site n and sites to the right of it (see Section 7.2 for details). It is determined 



by the configuration of the static disorder {cOn'} around the site, but not by the values of 



' This is a property of the two-dimensional diffusion. Indeed, as the chaotic region of the three- 
dimensional action space of the triple is concentrated around the straight line aii+gli = 0)2+ gh = <^3+gh, 
it is sufficient to determine the time it takes for the system to touch the corresponding region in the two- 
dimensional section, perpendicular to this line. Consider diffusion with the coefficient Z) on a square Lx L, 
divided into small cells a x a. It takes the time ~ a^/D to move to an adjacent cell, and the time ~ L?/D to 
cross the whole square. Thus, while crossing the square, the system visits ~ L^/cr cells, i. e., a fraction of 
the order of unity. In dimensionality higher than two, only a small fraction of cells would be visited. 
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the dynamical variables (over the latter, the thermal average is performed). It represents 
the probability that at a given instant of time the site n is hosting a chaotic spot, and 
plays the same role in the migration of a chaotic spot as the thermal weight e"^(^)/^ in 
the Brownian motion of a particle in an energy landscape ^(jc)]^ Let us estimate the 
typical value of w in the typical environment. Let Ng and A^^ be the number of oscillators 
involved in the guiding and layer resonances, respectively. Roughly, 



w~ y exp -^ exp --^-^ (T"pf^ (4.6) 

The first exponential is just the thermal activation probability g~(^*P^'o')/r for a resonance 



like that in Eq. (4.4), where the smallest frequency mismatch that can be obtained out 
of ~ e^« combinations of Ng frequencies (we do not distinguish between 2^«, e^«, or 
any other a^^ with a ~ 1), is ~ Ae"^s. The characteristic frequency of the separatrix 
motion is then Q. oc {TPp)^s^ where pNg is the typical order in r, needed to couple all the 
oscillators. The value of pis estimated in Sec. [7]as 1/2^/7^3. The second exponential 
corresponds to the exponential smallness of the Melnikov-Arnold integral e~^^^, where 
the smallest frequency mismatch of the layer resonance is cr ~ /S.e~^' . Finally, (rPp)^' is 
the strength of the layer resonance perturbation, determining the width of the stochastic 
layer. The leading exponential asymptotics of a sum of competing strong exponentials 
is determined by the largest term of the sum, or by the maximum of the exponent: 

1 (e~^^ e-^' Moll 
In - ~ min + +Ne\n U In^ In - , (4.7) 

W N,,Ni\ p {TPp)^s TPp) TPp p 

the optimal values being Nc ~ Ngln{l/TPp) ~ ln(l/p) ln(l /r'^p). Note that this already 
reproduces the main exponential part of the r, p-dependence of cr in Eq. ( |3.6a i. 



Thus, to each site is associated a random quantity w„. For different sites, n + n' , 
Wn and Wn' are assumed to be independenj^ Then, it is natural to apply the variable- 
range hopping picture, proposed for electron transport in strongly disordered solids by 
Mott ||52i : upon changing the guiding resonance, the chaotic spot typically shifts by a 
distance L which optimizes the compromise between the possibility to find a site with 
a large w,, (which is typically far away) and the exponential suppression of the hopping 
rate by the factor t^. It is also known that in one dimension the transport is determined 
not by typical hops, but by rare strong obstacles ll53]| . Such obstacles occur with a 



*We emphasize again that the sum 2]„ w„, equal to the average total number of chaotic spots on the chain, 
is dominated by rare "good" sites corresponding to resonant triples for t <c p, as discussed in Sec. |4.1[ Here 
we are interested in migration of chaotic spots between these triples. 

'This is not obvious a priori, given the extended nature of the guiding and layer resonances. Note, 
however, that replacing ci)„^ by (j>„|+i in Eq. 14.41 already changes the left-hand side by an amount ~ A. 
Thus, the optimal guiding resonances corresponding to sites n and n-l- 1 are likely to be completely different. 
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small probability, but in one dimension they block the transport efficiently as they cannot 
be bypassed, in contrast to higher dimensions. Following Ref. [611, we will use the 
term "breaks" to denote such obstaclesjj In the present problem, a break is a region 
of the chain where the chaotic fraction w assumes anomalously small values on many 
consecutive sites, i. e., a region of the chain rarely visited by chaotic spots. 

To implement these arguments, one has to know the distribution of the random quan- 
tities w„. The standard variable-range-hopping picture [52] assumes that ln(l/w) = Ais 
uniformly distributed. Here it is definitely not the case, so the estimate of the typical 
value ( |4.7| ) is not sufficient; one has to calculate the whole distribution function p{A). 
This implies calculating the statistics of different terms of the perturbation theory of an 
arbitrarily large order. The structure of the perturbation theory for Eq. ( |2.1| ) is discussed 
in Sec.[6j and p{A) is calculated in Sec. [7] The result is given by the Gumbel distribution: 

p^A) - -- exp f-Cipe[^'"'(i/^"^)]"'') . (4.8) 

dA \ I 

Here C is again a bounded function of Inp/lnr. It differs from C, which appears in 



Eq. (3.6a I for cr, by a factor of 2 at most. C\ is another logarithmic function of t and p; 



it enters Eq. ( 3.6a I in the combination ln( 1 /Cip), so it is beyond the logarithmic accuracy 



and has been omitted from Eq. (3.6ai. Eq. (4.8 1 is valid only for sufficiently large A » 



In (1/t''p), so it cannot be used to find the average value of w = e~ . Indeed, the average 
is determined by smaller A ~ ln(l /rp), corresponding to resonant triples, which should 
be treated separately. Given p{A) in the form (|4.8[), it is clear that the quantity which 



is uniformly distributed is the exponential under the derivative. Thus, AJ^ , the typical 
smallest value of A which can be found within a segment of a length L, is found as 

1 



L '-'""^ 



-Cip f^[Cin^(i/r"P)]-'^r - l)l , //P = Cln\l/TPp) In -^ . (4.9) 
\ /J ^ C]pL 



For L ~ I this agrees with the estimate (4.7 1. 



It is crucial that all the discussion of this subsection is valid regardless of the relation 
between r and p. Indeed, the typical order of the perturbation theory needed to form the 
guiding resonance is Ng ~ ln(l/p). Then the effective coupling between the oscillators 
of the guiding resonance, ~ (r''p)'"(^^^\ is always smaller than p, regardless of r. 

4.4. Macroscopic transport coefficients 

To relate the chaotic fraction w and the macroscopic conductivity cr, we use the fact 
that each segment between neighboring breaks contains many chaotic spots (this will be 
verified in the end of this subsection). Then, equilibration within each segment occurs 



** Other terms have also been used in the hterature for these objects, such as "blockades" 1641 or "weak 
Unks" |65l. The term "break" f61 66^ may not be the most precise, since the chain is not fully broken. Still, 
since the approach of the present work is analogous to that of Ref. j61J, the same terminology is used. 
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much faster than between different segments, so each segment k of the length Lk can be 
characterized by its own values of the temperature and the chemical potential, /i,t> T^, 
as well as the total action Itot,k and energy Hu- These latter ones slowly change due to 
the currents through the breaks J^, f^ (which can be viewed as the definition of the 
currents). It is convenient to denote the corresponding 2-columns by 



Qk = 



'tot,k 


, Jk = 


I'M 

I •'k \ 


, Vk = 


l^k/Tk 



(4.10) 



where the ^th break is assumed to separate the kth and the (k + l)st segments. Then 
energy and action conservation can be written as 



dQk 
dt 



Jk-\ - Jk- 



(4.11a) 



The current through the ^th break can be related to ^k, Tk,fik+i, Tk+i with the help of the 
diffusion equation in the space of actions {/„) (it is done in Sec.|9]), and can be written as 

-1, 



Jk-RViVk-Vk^i), 



(4.11b) 



to the linear order in the potential difference Vk - Vk+i- The validity of the linear response 
theory follows quite generally from the diffusive character of the system dynamics in 
the space of actions. This is equivalent to finiteness of the coefficient R^^, which is 
a symmetric 2x2 matrix and can be viewed as the conductance of the break. It is 
determined by the whole profile [A„} of the break, Rk = /?({/!„}), and is found in Sec. |9] 
Finally, we have a thermodynamic relation between Itot,k,Hk and Tk,Hk, discussed in 
Sec.[3]and |Appendix A [see Eqs. (A.la), (A.lbi], which we write here as 



Qk = LkQ{Vk). 



(4.11c) 



The explicit form of the function <3(V) follows from the discussion of Sec. [3] and Ap- 
pendix A} but is not important here. 



Now we pass to macroscopic equations, valid at large distances L » Lk- Eqs. (4. 1 la ) 
and (4. lie I lead to the continuity equation, dQ{V)/dt = -dJ/dx. To relate the current J 



to the voltage gradient, we note that in a stationary situation, when a constant uniform 



current J flows through the chain, the potentials can be found from Eq. ( 4. 1 lb i as V*: = 
- Yjk'<k Rk'J- The potential drop over a distance L including many segments determines 
the macroscopic gradient: 



dV 
dx 



V{x + L/2) - V{x - L/2) 



llkRk 



J = -o--^J. 



(4.12) 



L L 

Substituting the expression for the current in the continuity equation, we obtain the 
macroscopic transport equation: 



dQ{V) dV d ^,,^ dV 

= — cr{V) — . 

dV dt Bx dx 



(4.13) 




(b) 




X 




Figure 1: (a) A realization of A„ = ln(l/w„), where w„ is the chaotic fraction of the thermally weighted 
phase volume, associated with the site « (see details in the text, Sec. |4.3^ . Four breaks are shown, (b) The 
equivalent electric circuit with resistors and nonlinear capacitors (the bottom plate of each capacitor is 
grounded). 
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Eqs. (4.11ai-( [4.11c| ) are formally equivalent to those which describe an electric cir- 
cuit, shown in Fig.[TJb). Namely, each segment corresponds to a capacitor, each break to 
a resistor, while Qk and Vk correspond to electric charge and voltage. The only difference 
is that instead of the usual linear relation Qk = CkVt, where Ck is the capacitance, we 



have a nonlinear relation (4.1 Ici. In Eq. (4.13 1 dQ/dV plays the role of the macroscopic 



capacitance per unit length, and cr is the macroscopic conductivity. The physical reason 
for this analogy is that the theory of electric circuits describes nothing else but the trans- 
port of a conserved quantity (electric charge) driven by the gradient of the corresponding 
thermodynamically conjugate variable (voltage). In our problem the phenomenology is 
the same. 

Eq. (4.12 1 expresses the usual addition rule for resistors in series. Since cr^^ is 
proportional to the sum of resistances of individual breaks Rk, and these resistances are 
independent due to large spatial separation between breaks, then by the law of large 
numbers cr'^ is self-averaging at long distances. Thus, we can write 

CT~'= fR{{An})dP{{An}), (4.14) 

where dP{{An}) is the probability measure (per unit length) of the break configurations 



{An}, determined by product of independent measures (4.8 1. The integral in Eq. (4. 14 1 is 
mostly determined by the vicinity of a certain configuration which optimizes the strong 
competition between the exponentially growing R{{A„}) and quickly decreasing dP{{An}) 



(Fig. 15 1. /?({/!„}) is determined by the diffusion coefficient for oscillators in the break 
region. This diffusion coefficient contains the same small factors as the chaotic frac- 
tion Wn, SO they are strongly correlated. We take this into account by assu ming the dif- 
fusion coefficient to be proportional to e ^ ^ , 1 ^ C2 ^ 2 (see Sec. 8.4 1. Calculation 
of the integral in Eq. ( 4.14| ) is performed in Sec. |9} and is quite analogous to the opti- 
mization procedure of Ref. W\\ . The result is given by Eq. (3.6ai. The self-averaging 
of cr ' occurs at distances exceeding the typical distance between the optimal breaks, 
L* ~ e*"^" (1/t''p) Thus, it is L* that determines the border between the microscopic and 



macroscopic scales, Eq. (3.7 1. Note that this distance is indeed much greater than the 
typical distance between chaotic spots (a power of l/r, 1/p). 

4.5. Various remarks 

4.5.1. On the decay of correlations 

Let us see how the difference between oscillators which belong to a chaotic spot and 
the rest of oscillators manifests itself in their dynamics. One can study the two-time 
correlator: 

r'° dt' 

^(^t)= lim I — ^„(f' + 0^:(f'), (4.15) 



'a^°° Jo ^0 

or its Fourier transform, A„„((i>), which can be called the spectral function. For un- 
coupled oscillators (at t = 0) the correlator simply oscillates as A„„(?) = If,e~'^'^"'*'^'"^\ 
without any decay, so A„„(a>) is proportional to a 5-function. 
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For a chaotic spot residing on a resonant triple, the phase of each of its three oscilla- 
tors is randomized on the time scale of the destroyed separatrix, t^ ~ 1/Q, which should 
give the time scale of decay of correlations (we neglect the small residual correlations 
found in Ref. ll67ll for the standard map). Let us now consider an oscillator which is far 
away from a chaotic spot. Its phase is randomized due to fluctuations of the action /„, 
which translates into fluctuations of the frequency a>„(/„) = a>„ + gin- Namely, 

«A„(f' + t) ^nit') ex exp -/ r [oJn + ginit') + g dW)] dt" 

where 6ln{t") - In{t") - In{t') is the random increment of the action. For the diffusive 
dynamics 5In{t") ~ '^Dn\t" - t'\, where D„ is the diffusion coefficient. The correlation 
is suppressed when the random part of the phase becomes ~ 1, so it is natural to esti- 
mate the correlation time as t^ ~ {g^D„y^^^. The diffusion coefficient is exponentially 
suppressed at least as t^^, where L is the distance to the nearest chaotic spot, so the 
correlation time is exponentially longer than the time 1/Q for the resonant triple. 

However, it would be wrong to associate thus obtained tc with the decay time of 
A„„(f) and l/tc with the width of the spectral peak in A„„(w). Indeed, the peak is sup- 
posed to be centered at the oscillator frequency; but as the latter depends on /„, contri- 
butions from different times t' correspond to /„'s differing on the thermal scale T/\fj.\, so 
the peak is smeared over ~ gT/\fj.\ = pA. This smearing is at least of the same order as 
n ~ A ^Jmm{Tp,p^}. Thus, on the one hand, the spectrum of A„„(a») is broadened for all 
oscillators, and has no infinitely sharp 5-peaks, so the motion of all oscillators is chaotic. 
On the other hand, A„„((d) does not distinguish between oscillators belonging to chaotic 
spots and other oscillators. 

To eliminate this "inhomogeneous broadening" effect, one can use the conditional 
average 

/'" dt' Mt' + t) rnit') 6{In{t') - I) 

A„MI) - lim ^* . (4.16) 

j^° dt' dilnin - I) 

Then the main peak in A„„{a)\I) would be positioned at co = co,, + gl, and broadened 
by 1/n for oscillators belonging to resonant triples, or by an exponentially small value 
{g^Dn)^^^ for other oscillators. Note that the spectrum A„„(a»|/) also contains weak satel- 
lite peaks at different frequencies due to perturbative coupling to other oscillators. 

Another way to elimitate the "inhomogeneous broadening" effect is to consider a 
four-time correlator, 

r'° dt' 

Cnnnnit) = Um — lj/n{t' + t) ljjl{t') l!/*„(t') i/Znit' - t), 

to^oo Jo to 

where the regular oscillation simply cancels out. For uncoupled oscillators this quantity 
is simply given by I^, while in the coupled system it decays on a time scale tc which is 
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fast or slow for oscillators belonging or not to a chaotic spot, respectively, as discussed 
above. Accounting for the fact that a given oscillator can belong to a chaotic spot just 
for some intervals of time results in the correlator Cnnnnit) having both fast- and slow- 
decaying components with the corresponding weights. 

4.5.2. On the thermalization of oscillators 

Let us define thermalization of oscillators for a given trajectory of the system in 
the phase space {/„(?), <^n(0)- Let N oscillators n\, . . . ,nj^ he. called thermalized if the 
probability to find them in a given point of the phase space in the course of dynamics is 
given by the product of independent Gibbs distributions: 

(4.17) 
Here 6^ denotes the 5-function of a small but finite width s. 

Clearly, this must hold for oscillators which do not belong to a chaotic spot for most 
of the time, as their actions diffuse independently. However, this is not the case for os- 
cillators which actually provide the main contribution to the chaotic fraction. Consider, 
for example, a resonant triple on the sites n = 1, 2, 3 at r «; p. When it hosts a chaotic 
spot, the three actions stay close to the line oji + gh - a)2 + gh = 0*3 + gh, so their 
joint probability distribution does not split into a product. Even though the chaotic spot 
eventually leaves the triple, and the three oscillators can diffuse independently for some 
time, the relative contribution of such "non-Gibbsian" time intervals to the integral in 



Eq. (4.17 1 may be significant. 

Thus, the statement about thermalization of the chain is approximate. Namely, there 
is a small fraction of oscillators for which the Gibbs distribution is not valid. Never- 
theless, it does not affect the conclusions of Sec. |3] and Appendix A Indeed, when one 
calculates the total action and energy, their relative contribution is ~ min{Tp,p^), and 



thus is already beyond the precision of Eqs. (A.la), (A.lb). Still, the situation is some- 



what counter-intuitive: the "most chaotic" oscillators turn out to be the "least thermal". 

4.5.3. On the Kolmogorov-Arnold-Moser theorem 

Thermalization of the oscillator chain for general initial conditions, discussed above, 
seems to contradict the Kolmogorov-Arnold-Moser (KAM) theorem 1168116911701 . Ac- 
cording to the theorem, as the strength of the integrability-breaking perturbation r — > 0, 
the fraction of the chaotic part of the phase space must vanish. At sufficiently small r 
most of the phase space must be occupied by the invariant tori, where no thermalization 
is possible, while here thermalization is argued to be the case for general initial condi- 
tions. There is, however, an essential difference between systems with finite and infinite 
number of degrees of freedom. 
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As discussed in the previous subsections, the typical spatial density of the chaotic 
spots is M> ~ min{Tp,p^). This means that the probability for a chain of length L to 
contain no chaotic spots at all, is given by e'^^. This is just the fraction of the total 
phase volume, occupied by the KAM tori. Even though it goes to unity as t ^ or 
^ — > for any fixed L, if one lets L ^ oo first, it vanishes. Equivalently, even though 
one can prepare the initial conditions such that the chaotic dynamics is absent, they form 
a set of zero measure. 

Infinite number of degrees of freedom by itself is not sufficient to destroy the regular 
regions of the phase space. There are two more conditions: (i) disorder and (ii) exten- 
sive norm, and both are crucial. Indeed, early studies on one-dimensional systems of 
interacting particles have shown that even though the proportion of the regular regions 
tended to decrease as the number of particles A^ increased ItTI . it remained finite even at 
A/^ ^ oo as energy per particle was kept finite, provided that interactions were sufficiently 
short-ranged [72n73j|. At the same time, for disordered systems invariant tori have been 
shown to exist |l5l|74l[37], but they corresponded to a finite norm. 

One can view the collective effect of the infinite length, randomness, and extensive 
norm in the following way. Let us consider finite segments of A/^ neighboring oscillators: 
n = riQ + l,. . . ,nQ+N. Then, by simply varying hq one will always find a combination of 
frequencies a»„Q+i, . . . , aj„Q+A?, arbitrarily close to some given one, and the norm density 
on that segment will not be vanishingly small (i. e., bounded from below with a bound 
independent of «o)- 

4.5.4. On the effects of quantization 

The most essential element of the picture developed in this work is the chaotic spot - 
a collection of a. finite number of oscillators, whose frequency spectrum has a continuous 
component. Because of its the continuous spectrum, the chaotic spot can redistribute 
arbitrarily large amounts of energy between the oscillators (limited only by the total 
amount of the system's energy, which is infinite for the infinite chain), given enough 
time. In quantum mechanics, the energy spectrum of a system of any finite number of 
oscillators is always discrete, so the amount of energy which oscillators can exchange 
is finite even for infinite time. In such a system, irreversibility can arise only when 
the system has a possibility to visit an infinite number of discrete quantum states in 
an infinite sequence of tunneling events. A necessary (but not sufficient) condition for 
this is that the coupling between different states should be strong enough as compared 
to the energy mismatch between these states (this mismatch is always finite when the 
levels are discrete). The appropriate physical picture is that of Anderson localization- 
delocalization transition [IJ. However, now it should be considered not in the one- 
dimensional space of the chain, but in the space of excitations of the many-oscillator 
system. 

The quantum version of Arnold diffusion for a few oscillators was considered in 
Refs. f75ll76]| . The minimal number of oscillators giving the Arnold diffusion in each 
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model was taken, so the quantum system was effectively one-dimensional (moreover, 
of a finite length, since the number of oscillators was finite). Since in one-dimensional 
quantum problems with disorder all eigenstates are inevitably localized, the conclusion 
was that Arnold diffusion in a quantum system is always stopped by Anderson localiza- 
tion after a sufficiently long time. 

In the present problem the Arnold diffusion occurs in a system with many degrees 
of freedom, so one may expect the effective dimensionality of the quantum localization 
problem to be very high. In this situation, Anderson transition can take place. The 
control parameter for this transition would be the measure of "quantumness" of the sys- 
tem, i. e. the Plank constant h (of course, the measure of "quantumness" can also be 
expressed in terms of the physical parameters of the system at fixed h, as discussed be- 
low). Namely, for h not small enough thermalization and transport in the system can 
be completely blocked by Anderson localization. This transition would then be of the 
same nature as the many-body localization transition in a one-dimensional system of 
interacting bosons subject to disorder [77 1. 

Indeed, quantization of the Hamiltonian (|2.3|) corresponds to replacing the canonical 



variables 4/„, /(/'* by operators i/i'„, /i/',, with the standard commutation relation, [i/i-,,, /i^J] = 
ih. The resulting quantum Hamiltonian corresponds to the disordered Bose-Hubbard 
model. This model describes quantum bosonic particles, which interact with each other 
via contact potential, and whose motion is confined to an array of potential wells coupled 
by tunnelling. Then ^]j yh and i^„/ yh are creation and annihilation operators for the 
particle in the nth potential well, ha)„ is the energy of the bound state in the nth well, 
/itA is the tunnelling matrix element between bound states in the neighboring wells, 
and h?g corresponds to the Hubbard U, the interaction energy between two particles 
occupying the same potential well. The classical limit, formally corresponding to ^ ^ 
0, can be also expressed in terms of physical quantities. If one starts from the Bose- 
Hubbard model with a given disorder energy bandwidth W (corresponding to hA) at a 
temperature T, the classical limit at fixed W corresponds to T/W -^ oo, U/W -^ 0, 
UT/W^ = const « 1. 

The quantum bosonic system considered in Ref. llTTl corresponded to the continuum 
limit of the Bose-Hubbard model: the lattice spacing a — > 0, the site number \n\ -^ oo, 
maintaining x = na = const, and also t ^ oo, while keeping ra^ = const. Thus, no 
direct correspondence between the systems studied in present work and that of Ref. ITTTII 
can be established. 



4.5.5. On the Mott's law for variable-range hopping 



jtyp 



-d+CzW'/P 



Taking the typical value, A^ , from Eq. (4.9b, and maximizing T^e~^^ ^' l with 
respect to L gives the optimal distance for the hopping of the chaotic spot Lvrh ~ 
ln^(l/T^p)/ln(l/T), and the variable-range-hopping (VRH) estimate for o" in the spirit 
of Mott |[52l . This estimate is not supposed to be valid in one dimension because of 
rare breaks ||53]| . In our case, however, the only difference of this estimate from expres- 
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sion ( |3.6a| ) is in the factor C' , while the main functional form is the same (which still 
means that the resistance of a break is much greater than that of a long segment between 
breaks). This is due to the dramatic suppression (double exponential) of the probability 



p{A), Eq. (4.8 1, at very large A, which makes too strong breaks very improbable (in the 
standard Mott's picture A has a uniform distribution, corresponding to a constant density 
of states). 

It may also seem surprising that the picture of hopping transport, which usually 
gives the temperature dependence in the form of a stretched exponential {e~^^^ in the 
dimensionality d > I ||521 and e'^^^ in d = I ||53l ). translates here into e"'" '■^/'^\ The 
reason for this is that the hopping object, the chaotic spot, can be viewed as having 
many internal degrees of freedom. Indeed, on every site optimization is performed over 
a large number of available guiding resonances to minimize the frequency mismatch. 
As a result, the usual stretched exponential is stretched so much in this case, that its 
argument is logarithmic rather than power-law. This large number of internal degrees of 
freedom is also the reason for the strong suppression of p{A) at very large A. 

4.5.6. On the smallness of small parameters 

As mentioned in Sec.|3j the results of the numerical integration of Eq. ( |2.1| ) PTllQlfTSl 
[T6l |20l l2n |32l |33l |35l correspond to much faster dynamics than predicted in the present 
work. Moreover, the observed behavior could be understood assuming that upon ther- 
malization of the expanding cloud, the chaos uniformly spreads over all oscillators of the 
cloud [21 j. Again, this is quite different from the picture of rare chaotic spots. The av- 
erage density of the chaotic spots was estimated in the present work as vv ~ min|Tp,p^), 
and the numerical coefficient was not obtained. Thus, the chaotic spot picture should be 
valid when r and p are sufficiently small. The smallest values of r and p, for which the 
direct numerical integration of Eq. ( |2.1| ) produced reliable results on the packet spread- 
ing, are t = 1/40 andp about 1 [781 . So the question arises, starting from which values 
of T and p will the chaotic spot picture really work? 

One of the main physical assumptions used in the present work is that most of the 
normal modes of the linear problem are localized on one site of the lattice. Numerical 



studies of the linear part of Eq. ( 2. 1 1 have shown that at r = 1 the localization length is 
about 100, and it becomes ~ 1 only at r w 0.1 11791 [80l . This means that the disorder can 
be considered strong only for t < 0.1 at most. 

In Ref. [36], chains of finite length L were considered, and the probability for the 
chain to be in the regular KAM regime was analyzed numerically. At low densities, 
this probability was found to scale as exp[-10^(2p)^''^(2T)^^'''^(L/16)]. In the present 
work, it is interpreted as the probability to have no chaotic spots on the segment, e~^^, 
as discussed in Sec. 4.5.3 Thus, the results of Ref. fi^ approximately correspond to 
w X lO^T^p^. Apart from the difference from both w oc rp and w oc p^, obtained in 
the present work (the reasons for this discrepancy are unclear at the moment, so this 
issue deserves to be studied in the future), one should note the presence of the large 
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numerical coefficient 10^. Tliis means that for the chaotic spots to be rare very small 
values Tp < 10^^ are needed. 

These observations suggest that the transport regime, described in the present work, 
should set in at very small values of r and p, which probably have not been accessed in 
the numerical works up to the present date. 

5. Resonant triples 

This section represents a detailed analysis of chaos in a system of three oscillators 
whose frequencies happen to be close to each other. Although such resonant triples do 



not contribute to the main result for the conductivity, Eq. (3.6ai, they have the thickest 
stochastic layers. Estimation of their chaotic phase volume is the main task of the present 
section. 

First, we discuss a resonant pair of oscillators where chaos is generated by an ex- 
ternal perturbation, which can be studied by the standard methods ||54ll55ll59l . Next, 
the closed system of three oscillators is considered, for which no analytical solution is 
available. Still, an estimate for the chaotic phase volume can be obtained by analyzing 
various limiting cases and using the results for a perturbed resonant pair. 

The results of this section are relevant in the limit t <^ p. The opposite limiting case, 
T :i> p, requires analysis of high orders of perturbation theory, and will be discussed in 
SecO 



5.1. Two oscillators: separatrix 

The Hamiltonian of two oscillators has the form 



f/^-2/2.f'^ 



H = loih + Vi+ ^2/2 + f ^2 - 2tA V^cos(0i - 02). (5.1) 



Let us perform a canonical change of variables: 

<P = <f>i - (P2, l=^w^, 'P^^^^, l = h+h, (5.2a) 



^1+^-2,,^, 2,.,, .,.r,„K2 ...,^J^^2, 



H = / + 2. r + (wi - W2)/ + gl - 2tA V/V4 - /2 cos (p. (5.2b) 



The total action / is conserved, so only /, have a non- trivial dynamics: 

d^ dH ^ /cos0 

— - = — - - coi - (1)2 + 2gl + 2tA — , (5.3a) 

dt dl ^12/4 _ /2 

^ - - ^ = -2tA J/2/4 - P sin . (5.3b) 

dt d^ 

A separatrix in the phase space (/, 4>) appears when these dynamical equations have a 
hyperbolic stationary point. 
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Figure 2: A schematic view of the (/, ^) phase space, represented by the sphere. The circles = nil, 
^ i/> < 2n ("equator") and (p = 0,n, ^ & ^ n ("Greenwich meridian") are shown by dashed lines. 
The thick solid line represents the separatrix. The stationary points (one hyperbolic, at the separatrix self- 
crossing, and three elliptic) are shown by dark points. 



To find the stationary points, we make a change of variables / = (//2) c osQ, ^ 
& ^ n, thereby mapping the phase space onto a sphere (Fig. [2]). From Eq. (5.3b i it is 
seen that the stationary points require ^ = or tt. It is convenient to fix ^ = 0, and 
extend ^ < 2k, parametrizing the circle by only. Then we obtain an equation for 
stationary points: 



(jOl » » gl 

— - sin + T cos H sin 20 

2A 4A 



0. 



It has four solutions if and only if 

Otherwise, it has two solutions. To see this, we introduce 



r = 2. 



(wi - W2)2 + (2tA)2 



{giy 



00 = arctan 



2tA ' 



(5.4) 



(5.5) 



(5.6) 



so Eq. ( 5.4 1 becomes sin 20 + r cos(0 - 0o) = 0. At the borderline, it has three solutions, 
one of which is degenerate. Equating to zero sin 20 + rcos(0 - 0o) together with its 
derivative, we obtain 



r^ = 1 + 3 cos^ 20 = 



1 + tan^ 



tan(0-0o) = -2cot20. 



(l+tan2 0)3' 
From the second equation we obtain tan© = (cot0o)'^^. Substituting it into the first 



equation, we obtain Eq. (5.5 1 
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Thus, for the separatrix to exist, it is absolutely necessary to have 2tA < gl. This is 
the reason for making the assumption t «; p, as the typical value of / ~ gT l\iJ.\. Given 
tA «: gl, the condition (5.5 1 can be taken in zero approximation, \oj\ - tL»2l < gl- This 
is si mply the condition for existence of /i,/2 > satisfying co\ + gl\ = W2 + gh, see 
Eq. (4.2 1. Also, at r <s g/ the separatrix is concentrated in a narrow region of 7 around 
the resonant value I - (a»2 - a>i)/(2g), so in the last term of Hamiltonian (5.2bi one can 
simply set I = (0*2 - a>i)/{2g). Then the Hamiltonian becomes equivalent to that of a 
pendulum, 



//- 



9 / \ 2 9 

a>i + a»2 g 2 i^i ~ ^2) , ^If C02-C0]_\ Q sgnr 



■I+-r 



+ )i I 



2g 



2 4 4g 

and the pendulum frequency is given by 

5.2. Two oscillators: stochastic layer 



2g 



COS( 



(5.7) 



(5.8) 



Suppose now that a time-dependent perturbation is added to the Hamiltonian (5.1 1 



H' = 2V{I\, I2) cos(m\4>i + m2(p2 - ojt) 



(5.9) 



The frequency cu is a linear combination of frequencies of other oscillators. Their motion 
is taken in zero approximation, so effectively they exert an external time-dependent force 
on the resonant pair. This force produces a stochastic layer around the separatrix, as is 
discussed in detail in Refs. 



Perturbation (5.9 1 leads to a change in the old integrals of motion, I and H, 



dl 
dt ~ 
dH 
dt 



dH' 
dHdH' 



dH dH' dH dH' 

+ 



dl d<p dl d(p d(p dl 



(5.10a) 
(5.10b) 



Out of the three terms on the right-hand side of Eq. (5.10bi, the first one is essentially 
the repetition of Eq. (5.10a), while in the third term dH/d^ = 0{t), so it can be ne- 



glected. Being interested in the perturbation of the pendulum motion in the vicinity of 
the separatrix, one can take the values of l\ , I2 at the point of resonance, the unperturbed 
dependence of (^ = {dH/dI)t, and 0(f) right on the separatrix: 



0(0 ^ ± [4 arctan e^^'~'°^ - n] . 
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(5.11a) 

(5.11b) 
(5.11c) 



The two signs of ^{t) correspond to the two branches of the separatrix (Fig. |2]l. The 
changes in / and H during one passage of the separatrix are given by 



61 = 



CO 



2{mi + m2)V I sin 
2V 



nil — 1712 



4>{t) - AQt 



dt = 



6H = 



Q. 

dH 



{nil + '«2)'5^±(m,-m2)(A) sin(Aafo), 



(5.12a) 



dl 

2Vco 



SI + (mi - m2)V 



J 



4(0 

dt 



sin 



m\ — m2 



4>{t) - Ant 



dt = 



a 



■ ^±(mi-m2)(A) sin(AQfo), 



(5.12b) 



where the Melnikov-Arnold integral Jl^iA) is defined as 



^m{A) 



CO 

J 



exp 



im 



-(4arctane* - n) 



e-'^'ds, 



(5.13) 



A 



CO - (mi + m2)(a»i + C02 + gI)/2 



Since J?l_„,(A) - ^m(-A), and J?l,„>o(A < 0) - {-iy"Jl^{\A\)e-"^^\ for |A| » 1 and for 
a given perturbation (5.9 1, one of the branches of the separatrix is much more efficient 
than the other one. If we fix t > 0, mi - m2 > 0, and A > for definiteness, the main 
contribution comes from the upper branch of the separatrix, shown in Fig. [3] At A » m 
the Melnikov-Arnold integral is exponentially suppressed, 



\m-l 



C2A)" 
:?1™>0(A » m) « 47r ^^ e-'''"\ 
r{m) 



(5.14) 



where r(m) is the Euler's gamma function. 

Let us represent the pendulum Hamiltonian [the last two terms of Eq. (5.7i] as 
(Q?/2g){h + 1), where h is the dimensionless energy of the pendulum, counted from the 
separatrix. Then h measures the distance of a trajectory from the separatrix. Denoting 

(5.15) 



h, = ^ A^ ^„„_^2(A), 9 = AQto, 



Q2 



and using Eqs. ( 5. 12a 1, ( 5. 12b 1, we write the change in h upon passing the upper branch 
of the separatrix as 

h \-^ h -\ sin 9, (5.16a) 

A 

For trajectories close to separatrix, |/i| «; 1, the time spent in the vicinity of the hyper- 
bolic point diverges as {l/Q.)ln{32/\h\). This gives the value of the phase 9 for the next 
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|- 0)2-0)1 




Figure 3: Stochastic layer around the separatrix in the (7,0) phase space, as shown by the shaded area 
(t > 0, mi - m2 > 0, and A > are assumed). 



passage of the separatrix: 



32 
+ (2)Aln-. 
\h\ 



(5.16b) 



The factor of 2 should be taken for trajectories inside the separatrix (pendulum oscilla- 
tions), since the hyperbolic point is passed twice before the system returns to the upper 
branch. For trajectories just above the upper branch (pendulum rotations) there is no fac- 
tor of 2. Trajectories below the lower branch are excluded from consideration since they 



never get to the upper branch. Eqs. (5.16ai, (5.16bl define a mapping which is often 
called separatrix or whisker mapping. If the mapping is iterated starting with a suffi- 
ciently small value of h, the phase 6 tends to randomize. For |A » 1| the mapping can 
be linearized in h around one of the resonant values /i„ by introducing a new variable z: 



h — hn 

z = (2)A — ; , hn - 32 exp 

h„ 



Inn 
(2)A 



(5.17) 



The separatrix mapping then reduces to the so-called standard mapping, also called 
Chirikov-Taylor mapping: 



z^-^ z + KsinO, 6 t-^ 6 - z, 



(5.18) 



where the stochasticity parameter is defined as K{h„) - {2)hjh„. As discussed in 
Refs. ll54l l55l [5 9 J . the standard mapping becomes chaotic when K > Kc. Numeri- 
cally, Kc is very close to 1, and we neglect Kc - 1 in what follows. This determines the 



stochastic region for the separatrix mapping (5.16aM5.16bi as \h\ < {2)hs. Using this 

29 



result, one can calculate the total phase volume of the stochastic layer (the shaded area 
inFig.[3]l: 

J 2n Ing \ ' 2hs 2 h^ j 

layer 



w. 



For |A| » 1 the change in h upon one iteration is much smaller than the layer width, 
so the dynamics of h is diffusive. If the phases at successive iterations are assumed 
uncorrelated, the average square of displacement in h after A^ iterations is {N/2)(hJA)^. 
It takes the time N6t, where 6t = {2)Q.~^ ln(32/(2)/j) is the time step, which gives the 
following estimate for the diffusion coefficient in h: 



Dm = - ^ 



(2) In 
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(5.20) 



In fact, the assumption of completely random phases is never fully correct. For the stan- 
dard mapping, Eq. ( 5. 1 8 1, even at relatively large K ~ 10-20 small residual correlations 



of the phase (at the level of 10"^ - 10~^) have been found f^l\. These residual correla- 
tions are usually explained by the system sticking to various structures inside the chaotic 
region, such as stabihty islands where the diffusion is anomalously slow, or accelerator 
modes, where the diffusion is anomalously fast. Still, in spite of the residual correla- 
tions, the diffusion coefficient in z for the standard mapping for A" > 2 was found in 
Ref. 1671 to agree quite well with its random-phase value K^ 12 together with the lead- 
ing in 1/ V^ correction lISTl : the largest relative deviations from K^ 12 were less than 
100% and occurred for special values of K corresponding to accelerator modes. The 
accelerator modes should be much less important for the separatrix mapping than for the 
standard mapping; indeed, since K - K{h), the fast diffusion in h will quickly move the 
system away from the special values of K. 

The random phase assumption is completely wrong close to the stochastic layer 
boundary. Indeed, the system cannot leave the stochastic layer and end up on a sta- 
ble periodic orbit, so the diffusion coefficient must vanish at the boundary. According to 
numerics, it vanishes as a power law, D/,(/j) oc ((2)/j.v - \h\f-^^- |[54| . As will be seen later 



(see Eq. ( 8.28a i and Appendix Fl, what will matter for the present problem is a certain 

average of the diffusion coefficient over the stochastic layer, so we estimate the relative 

error of Eq. (5.20 1 to be ~ 1 , which is beyond the accuracy of our final result, Eq. ( |3.6a I. 

As the perturbation causes changes in /, the latter also undergoes diffusion. However, 



Eqs. (5.12ai, (5.12bl, rigidly relate the change in / to the change in h on each step: 



(mi-l-m2)n 
61 = -— 5h, 

2gA 



(5.21) 



so, as h is limited to a thin stochastic layer, a single perturbation term (5.9 1 cannot 
produce a large excursion in /. 
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Figure 4: The allowed region of /i,/2, shown by the shaded area. The resonance manifolds defined by 
Eqs. l|5.24a^, l|5.24b|l are shown by dashed lines. The ellipse represents a constant-energy surface. 



5.3. Three oscillators: resonances 

The Hamiltonian of the triple is given by 

3 2 

7/ = ^ Y^Jn + 2 ^" ~ ^"^^Xi V^«^n+1 COS(^„ - 0„+i). 
n=\ n=\ 

We perform the following canonical change of variables: 

(p\ + (f>2 + <p3 
<P - ^ > (p\ = (p\ - (f>2 , 02 = 02 - <^3 , 

3 3 

a)\+a)2 + oj-i . gf- 

tl — I -\ h 

3 6 

+ (Wi - 0J2)h +gi\ + (W2 - ^3)72 + gll - ghh - 



(5.22) 



(5.23a) 
(5.23b) 



- 2tA ^Jhl2 cos (p\ - 2tA ^jl2h cos 02 • 



(5.23c) 



The allowed region of /i,/2, determined by the conditions I\,l2,h > 0, is a triangle 
-7/3 <Ii < 27/3, 1 1 - 7/3 < 72 < 7/3, shown in Fig.li 
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We consider two possibilities for the guiding resonance: either a»i + gh =0*2+ gh 
or a»2 + gh = <^3 + gh- The resonance between the oscillators 1 and 3 has the effective 
coupling ~ T^, so its stochastic layer is negligibly thin. The resonance conditions are 



^^0 ^ h-2h = '^^^^^ 
at g 

at s 



They intersect at the point 



y. —2(jj\ + (i>2 + a>3 y —a)\ - (x»2 + 2(i>3 

^Ir - T , hr = — 



3^ 



3^ 



(5.24a) 
(5.24b) 

(5.25) 



This point lies inside the allowed triangle if gl > 2a>j - coj - a>k for any {/, j, k] being a 
permutation of {1, 2, 3}, or, equivalently, 



gl >3 max {(On) - V W„ = glmin- 



(5.26) 



n=\ 



The importance of the intersection point can be seen by considering the result of the 
previous subsection for the 1-2 resonance, where the separatrix-destroying perturbation 
comes from the coupling to the third oscillator. Let us set mi = 0, m2 = 1, w = W3 + gh, 
V = -tA ^/hh in Eq. (5.9 1, and take the large-A asymptotics (5.14i of the Melnikov- 



Arnold integral appearing in Eq. (5.151 for h^, which, in turn, enters Eq. ( 5.19| ). This 
gives the following estimate for the 73 -dependent volume of the stochastic layer around 
the 1-2 separatrix (numerical and logarithmic factors are omitted): 



w: 



(12) 



A 

a 

/(12) 



lhhTK^e-''''l\ 



A 



\oj2 -UJ3+ gjh - h)\ 



(5.27) 



If we now integrate W^ over 73 to find the total volume of the stochastic region in 
the space ih,h, h), the exponential factor will ensure that the integral is dominated by 
values of 73 close to the resonance intersection point. 



5.4. Three oscillators: chaotic phase space 

Let us expand the Hamiltonian, Eq. ( |5.23c i, around the resonance intersection point, 
Eq. (5.25 1, writing 7] 2 - ^[2 "'' P^-^' 



H = 77,(7) + 77(pi,p2;0i, 02) = 

(a>i - 0*2)^ -I- (a>2 - (x»3)^ + (a»i - 61^3)^ wi -1- 0*2 + ^3 , gl^ 
6^ ^ 3 ^"6" 

+ gpl - gPiPi + gp\ - 2Vi cos 01 - 2^2 cos 02 ■ 
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(5.28) 



Taking Vi^ at the intersection point is allowed when the resonance intersection point 
(/[./p is not too close to the sides of the allowed triangle in the (/i,/2) plane. This is 
possible when gl » tA. 

The main task of the present section is to estimate the fraction of the total (thermally 
weighted) phase volume, occupied by the chaotic part: 



/ 



-{H-tiI)/T 



n„=i dlndcPnliln) 



w = 



chaotic 



je-(H-Ml)/T Yll^^dlndcPnliln) 



(5.29) 



At gT l\^i\ «; A «; \^\, the denominator can be simply taken to be (r/|//|)^. Inside the 
allowed triangle the typical value of H ~ gP «; T, so its contribution to the thermal 
weight can be neglected. Then the calculation of w can be done in two steps: 



Wssil) 



dl\d^i dl2d^2 

In In 



(5.30a) 



w 



-! 

chaotic 

3 ^ 

= 1^1 Ce-^"'^'^->^'^'''WAI)dI. (5.30b) 



Not being able to calculate Wss(I) at arbitrary Vi, V2, we first consider the case when 
one of the couplings is much stronger than the other (e. g., Vi » V2), and then extrapo- 
late to Vi ~ V2. At V2 = we write H as 

H\y^^o = s(p'-y) ~ ^^' ''°' h + l8Pl^Hi + ^- gpl , (5.31) 



and find the separatrix solution for ^\, as in Sec. 5.2 



0i(O = ± [4 arctan £"'<'-'») - n] , hit) = -'^ + ^- gpjt, a, = Vw". (5.32) 



Again, we calculate the changes in the action p2 and the pendulum energy H\ : 

V2 



Sp2 - - I 2V'2 sin 02(0 dt = 



IWl 



■ ^±\ 



^gP2 \ . 'igPltQ 

Ism 



A^|gV\l 4V^ 



dHi = 



-! 



—- IV2 sin 02(0 dt = -11^ j{^^ 

OP2 



2V^ 



^gP2 \ . ^gPlto 

Ism 



^VgVj 4^/gVi 



(5.33a) 
(5.33b) 



where ^±i(A) = ±27re*''^''2/ sinh(;rA). 

Again, one arrives at the separatrix mapping, analogous to Eqs. (5.16al-(5.16b 1. The 
main contribution to the stochastic phase volume comes from g\p2\ ~ ^gV\, i. e. |A| ~ 1. 
This is still much larger than the inverse time of the passage of the separatrix due to the 
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logarithmic factor, so the phase is randomized after a few iterations. However, at |A| ~ 
1 the change in the Hamiltonian per iteration is not small compared to the stochastic 
layer width, so the standard approach for determining the layer width (linearization of 



Eq. (5. 16a I around some value of h and reduction to the standard mapping) does not 
apply in this case. The diffusive picture for the motion across the stichastic layer is not 
valid either. Also, for the pendulum rotation, H^ > IVx, the contribution oc g'^^l'^l/^ 
from the "inefficient" rotation direction should also be kept. Thus, the stochastic layer is 
determined by the inequalities 



■ AV^K^^Ro^t^) <Hi- IVi < 2y2A^^±(A), 



(5.34) 



where J?lo,±(A) are some functions, analogous to Melnikov-Amold integral (5.13 1, about 
which we only know that 



J{o{-K) = J?lo(A) , :?l±(-A) = :n^(A) , (5.35a) 

|A| ~ 1 : :no,±{^) ~ 1. (5.35b) 

|A| » 1 : :noiA) w Ine'''^^^'^, :5^±(A) ~ Ine''^^^-''^'^. (5.35c) 



Now we can calculate the volume of the stochastic layer in the pendulum phase space 



analogously to Eq. (5. 19 1 for each p2, and then integrate over p2 [equivalently, over A 

3m/(4ViW)]: 



Wv 



37T g 



! 



dA 



4A^:^.(A)lnl^ri^ . y A^:^.(A)ln^^ri/^ 
A2^„(A) 4" A2j?l^(A) 



(5.36) 



The integral is dominated by |A| ~ 1 where the functions ^o,±(A) are not known. Thus, 
one can only conclude that 

Vi Vi 
W,, = C, — In -^ , (5.37) 

g V2 

with some unknown constant C^s ~ 1 ■ 

At V2 ^ Vi, the main region of integration over p2, \p2\ ~ y/VyJg, is far away from 
the separatrix of the {p2, ^2) motion, which occurs at \p2\ ~ ^JV2/g■ When V2 ~ Vi, the 
two separatrices should be treated on equal footing, and the author is not aware of any 
accurate method for description of the stochastic layer. A reasonable extrapolation of 



Eq. (5.37) is 



''ss ~ ^ss 



Vrr, 



In 



C V. 



ss ' max 



(5.38) 



where V,nax = max{Vi, V2), V,„,„ = min{Vi, V2), and Css,C'g^ ~ 1 are some bounded 
functions of Vmax/Vmin- It is worth reminding that this expression is valid when the reso- 
nance intersection point is inside the allowed triangle and not too close to its boundaries. 



i. e., / satisfies inequality (5.26 1 with a positive correction ~ rA to its right-hand side. 
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Now let us turn to the integral in Eq. (5.30b I. The whole / integration range consists 



of three parts: (i) the smallest I, < I < I ■ , when the allowed triangle is so small 

m . . -Mine) 

min 



that none of the dashed lines in Fig. 4 enters it, (ii) intermediate /, /^.^ < I < I,„in 
when one or both dashed lines enter the triangle, but not their intersection point, and 
(iii) sufficiently large / > /„,,„, when the triangle contains the resonance intersection 



point [Eq. ( |5.26[ )]. For a typical configuration of the random oscillator frequencies, 

A/g. Hence, even though Wss in the intermediate region has 



Mine) 



'^mirit ^mm 



Mine) j j 

^min ^^min,h 

an additional smallness ~ e^^^^ from the Melnikov-Arnold integral, a thermal factor 
~ e^^^ is gained with respect to / > pP"'"'' due to the higher activation barrier of the 
latter. However, as will be seen below, the main contribution to the disorder average 
of w comes not from the typical configurations, but from those where /,„,„ ~ T l\ij.\. In 
this case the contribution from the region / < I < Imin can be neglected, and one can 



focus on / satisfying inequality (5.26 1 and use Eq. (5.381 for W., 



As 1 increases, the intersection point enters the triangle through one of its sides, 
depending on which of (x>i, a>2, a>3 is the largest. It is convenient to rewrite Hr{I) as 



HriJ) = Ea+ max {U)„ -//)(/- Imin) + 7 (^ " ^min) , 
l^Hs£3 6 



Ea = 3 max 



(w„ - fl) 



l^n^3 2g 



2 A r... ^a2 
2g 



2 

n=l 



{co„ - IJ.f 



(5.39a) 
(5.39b) 



The contribution of the last term in Eq. ( 5.39a i to the thermal weight can be neglected. 
Ea can be viewed as the activation barrier for the triple. Its probability distribution is 
calculated as (using A «c |//|) 



Pcj{Ea) 



I/^|A 



A/2 

/doji 



dco2 dco3 



A A 



-A/2 



3 max 

Kns£3 2g 



-z 



(Ojn - H) 



0, 



n=\ 



2g 



-Ea 



'^(l-^^), 0<^<1, 
1-l^f, 1<^<2, 

l\^i\A) MA 



(5.40) 



otherwise. 



Let us assume that wi > W2,W3, which occurs with probability 1/3. Denoting g{I - 
Imin) = oj', 0J\-iO2 = cou, oj \ - LOj, - a>i3, wc obtain 



ICsstA _£ 
w - e ' 

'i{gT|\^l\)^ 



'"1'^ ^ ^lu)\oj' + 3con)e-''"^''^«'^hn 



C'^^ VW + 3aJi3 



dio'. (5.41) 
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The integral can be calculated analytically in two limiting cases: 

InC;^, wi2,wi3 «;gr/|//|, 



w 



3 gT ^ 



' In '■' _ , aJn,0Ji3»gT/\^i\, (5.42a) 



16 gT 



gT 



where y = 0.577 ... is the Euler-Mascheroni constant. For cot, > a)i,a)2 one obtains the 
same result with a>\ and a>j interchanged. For a>2 > wi > oj^ the calculation is fully 
analogous, and one obtains 



w 



2C,,T |//|A _j;jj 

3 gT "" 



( InC',, 



W21,W23 <^gT/\lj\, 



37r a;2^2i , (C;.)'^23 ^,, , (5.42b) 

16 ^r (JJ21 



Again, for W2 > <x»3 > a>i it is sufficient to interchange coi and W3. 

Now one can average w over the disorder, which gives the average spatial density of 
the chaotic spots: 



w 



A/2 

/ 

-A/2 



doji da>2 da>3 

'a a A~ 



W = C„ T 



l/^|A 



CwTp , 



(5.43) 



where C^ is a constant. The average is dominated by rare triples which have Ea ~ T, 
that is |aJi2| ~ IW23I ~ gT/\iu\. Since expressions (|5.42a[)-([5.42b|) cover only limiting 



cases, no analytical expression for the constant C„ is available. 

The estimate ( 5.43 1 is not very much sensitive to the fact that in Eq. ( |5.29 1 the chaotic 
fraction is defined using the thermal weight. The latter only fixes the typical value of the 
action /„ ~ gT/\ij\. Let us repeat the estimate simply fixing h - h - h - ll'i, so instead 
of the thermal exponential in Eq. (|5.29|) we put 6{l\ - 1/3) 6{l2 - 1/3) dih - 1/3). Then 



V\ -'V2- (2/3)tA/, and w is simply given by the probability that the point Ii = I2 = 
hes inside the area W,,v(/) ^ C^^rM/g (we denote C"/ = (2/3)C,, In C^,) around the 
resonance point (5.251, which gives vv - 3C'//Tg//A. Defining p = ^//(3A), we obtain 



an estimate fully analogous to Eq. (5.43 1. 



6. Perturbation theory 

6.1. General remarks 

One could formulate the peturbation theory for Eq. ( |2.1| ) by taking independen t os- 
cillations, i/',, (f) = -^fj^ e''^^"'^^^"'"'''"^" , as the zero approximation and solving Eq. (2.1 1 



for i/^nit) by iterations. Clearly, such perturbation theory is divergent. Indeed, if it were 
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convergent, there would be no chaos. Formally, its divergent character can be seen 



already from the exact solution for two oscillators, given in Sec. 5. 1 The pendulum fre 



quency depends on t as Q oc Y|t| [Eq. (5.8 1], while the dependence on g has a threshold 



[Eq. (5.5 1]. As the dependence on t and g is non-analytic at T,g -^ 0, the perturbation 
series has zero radius of convergence. 

The crucial point of the present work is that instead of studying the whole perturba- 
tion series, the most dangerous terms of this series should be identified (more, presicely, 
the probability for them to occur), and the subsequent analysis of their effect should 
be performed non-perturbatively, using an effective model, such as driven pendulum. 
In other words, once a small frequency difference in the denominator has occurred in 
the pertubation series, one step back should be taken, the combination of the oscillator 
phases corresponding to this frequency combination should be called the phase of the 
effective pendulum, while the prefactor multiplying the inverse of this frequency dif- 
ference should be called the effective coupling, which multiplies cos in the effective 
pendulum Hamiltonian. 

Divergence of the perturbation theory due to resonances corresponds to an instability 
at long times. Indeed, perturbation theory in the time domain in the nonlinearity has 
been shown to converge only for not too long times 10711221 . Its divergence at long times 
signals the chaotic instability. In fact, the perturbation theory developed here is closely 
related to that of Refs. |[T7ll22l . as will be discussed in more detail the end of Sec. |6.3| 

The perturbation theory developed here is designed to generate the effective many- 
oscillator coupling in high orders in r and g. We define the effective coupling between 
2N oscillators ni,. . . ,nN,ni, . . . ,hi^ (these are, generally speaking, 2N arbitrary inte- 
gers) as the coupling term in the Hamiltonian which involves no other oscillators and 
leads to the same equations of motion for these 2N oscillators, as the equations obtained 



by eliminating all other oscillators from the full equations of motion, Eq. (2.1 1. Thus, 
the effective coupling is the sum of terms of the kind 

y'"' = ^<?!L.,n.....'A:. ■■■C'A.: ■■■<A.„, (6.i) 

Strictly speaking, the above definition of the effective coupling is meaningful only ex- 



actly at resonance, a*,,, -i- . . . -i- a»„_^ = cofn -i- . . . -i- a>jj^, as will be seen in Sec. 6.3 This 
does not pose a problem, since it is just in the vicinity of the resonance that this coupling 
is needed, where it determines the parameters of the effective pendulum. Thus, the fre- 
quency difference laj,^ -i- . . . -i- o),,^ - a*,!, - ... - oj^J is assumed to be smaller than any of 
the intermediate frequency denominators involved in the coefficient 7C. In other words, 
the 2A^-oscillator perturbation term is generated under the assumption of no resonances 
for 2N' < 2N oscillators, involved in the construction of this 2A'^-oscillator term. Indeed, 
if a 2A'^'-oscillator resonance occurs, it would be more efficient in producing chaos due 
to larger effective coupling obtained in a lower order of the perturbation theory, and then 
V*-^^ is not needed at all. 
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Since we are interested in the leading approximation to the effective pendulum fre- 
quency, we will take the lowest order in g which couples the 2N oscillators altogether 
(that is, g^~^). Similarly, we will take the lowest possible order in t, which is determined 
by the spatial arrangement of the oscillators. According to this argument, one should ig- 



nore terms of the kind ( |6.1[ ) where some of n's coincide with some of h's. Indeed, the 
coinciding ones do not enter the slow phase 0„, -i- . . . -i- 0,,^ - 4>ni - ■ ■ ■ - <PnN^ ^o the 
same effective pendulum can be constructed in a lower order of the perturbation theory. 
Thereby we neglect all secular terms of the perturbation theory, and, consequently, the 
nonlinear frequency shifts in the frequency denominators involved in the coefficient 7C. 
Due to the condition g\iJ/n\^/A « 1, Eq. (2.121, these shifts are small compared to the 



dispersion of the random bare frequencies a>„, and thus give small corrections to the 
statistics of the denominators. Such corrections are beyond the accuracy of the present 
work. 

6.2. Diagrammatic rules 

In this subsection we describe the formal procedure for generating 2A'^-oscillator 



couplings (6.1 1. Its derivation will be given in Sec. 6.3 

For two given oscillators n\,n2, the coupling by tunneling is taken in the lowest order 
in T (that is, tI"'""^!). This is done by defining the Green's function (propagator) as 

Gn.mioj) = i (-tA) i (-tA) . . . (-tA) ^ , (6.2a) 

(jJ — COn ■ (jO — Cl)„ +1 (Ji) — (jJyi 

where «mm = min{«i, «2K '^m«.r - max{?ii, ?i2}- We also define the shortened propagators: 

G'n^n2i(^) ^ (^ - w,,,) ^^[^^(w), G," ^^(oj) ^ Gn^mioj) {oJ - w,,,). (6.2b) 

This is nothing but the Green's function of the linear problem, taken in the leading order 
of the locator expansion. In contrast to other works where the solutions of the linear 
problem are assumed to be known and their properties are postulated |[T7ll20ll2Tll22l|23l 
[3T]|33l|37l, here these solutions are constructed explicitly using the perturbation theory 
in T. This is the reason why the Green's function of the linear problem appears. 
The effective 2A'^-oscillator coupling is obtained from the following procedure. 

1 . Draw all distinct diagrams for the left-going solid line obtained after N - I itera- 
tions of the diagrammatic equation 



(6.3a) 



supplemented by an analogous equation for the right-going line 

A-, — = 



-^ ' . (6.3b) 
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To illustrate the diagram generation, all distinct diagrams for A^ = 2, 3, 4 are shown 
explicitly in Fig. [5] 

Each of the free ends of solid lines corresponds to some of the sites «i, . . ., n^, 
hi,. . . ,nN. Each dotted-solid-solid vertex also has a site index, the two ends of 
the same dotted line corresponding to the same site. 

Each solid line has a frequency argument. For external lines it is the frequency of 
the site at its free end. For internal lines it is determined from the conservation 
law holding for each dotted line: the sum of the two frequencies of the entering 
solid lines is equal to the sum of the two outgoing frequencies. An example of 
labeling a diagram is shown in Fig. [6] Note that this procedure is consistent only 
when aj„^+ ... + aj„^ = ajni+ ... + o^hn- 

One obtains an analytical expression corresponding to each diagram by (i) associ- 
ating a factor g to each dotted line; (ii) associating a factor Gnn'{(^) to each internal 
solid line which starts on site n, ends on site n' , and carries frequency co; (iii) asso- 
ciating a factor G'„^„{a>ni,) or G'4.(<x>,ij.) to each external outgoing or incoming line, 
respectively, with free end on site «/t or «jt; (iv) summing over internal site indices. 
For example, the diagram in Fig. |6] produces an expression 



5. Summation over all diagrams gives the coefficient 7C„ „ ^ ^ in Eq. (6.1 1. 



6.3. Derivation of the diagrammatic rules 



Let us write the equations of motion, Eq. (2.1 1 in the form 



{id, - aJn)ll'n = -TAilf/n+l + ij^n-l) + ^n 2_j Snn'l^l'll'n' , (6.4) 

n' 

where gnn' - gSnn' is introduced to simplify the diagrammatic rules (otherwise, ad- 
ditional combinatorial coefficients would have to be assigned to each diagram). The 



right-hand side of Eq. ( 6.4 1 can be viewed as an external force acting on the oscillator n. 



Then, the response of the oscillator to an external periodic force is determined as 

(idt-conji/yn-yjae-'''"' ^ ^„ = ^f e"'""' + V ^^^^ . (6.5) 



a 



Let us start with two oscillators, n and n (assume n < nfoi concreteness). Suppose 
the «th oscillator performs some given motion 

-^M(^)e''''\ (6.6) 
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N=2: 



N=3: 
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N=4: 
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+ 
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..< 



/<CJ 
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Figure 5: All distinct diagrams generating effective couplings for A' = 2, 3, 4 oscillators. 



n, 



K 




CDn + Ww - W^ 



-n, 



CD^/ CD^^+ CD^^- CD„^- CD„^ 



CD; 



Figure 6: Labeling the lines and vertices of a diagram. 
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and let us find the corresponding force acting on oscillator n. To the leading order in 
T (which is n - n) and g (which is zero), this force is found by writing Eq. (6.4 1 in the 
Fourier space, and eliminating oscillators n + I, .. . ,n - I perturbatively: 



-^ (-tA) (-tA) . . . (-tA) (-tA) ipnico) e-''^' ^ 

In 0) - COn+i OJ-OJn-l ,^ „. 

In 



r 



We are interested in the case when the spectrum of ij/fj is concentrated around co v a»,-, w 
ojn, i. e. slow modulations of the periodic oscillator motion. Then one can replace 
G'^'fjico) -^ G"^{(jOn) and take it out of the integral. Then the effective equation of motion 
for the «th oscillator is 

(id, - OJn) ifJnit) = G'^iaJn) ifJh{t) , (6.8) 

which can be obtained from the effective coupling Hamiltonian 

V = G'2{a)-n) ipl^ph- (6.9) 

This corresponds to the diagram consisting of a single solid line with two free ends, so 
we have introduced the Green's function G'" , shortened on both ends. To the coupling 
Hamiltonian in Eq. (6.9 1 another term should be added, which has the same form, but 
with n and h interchanged. The resulting Hamiltonian is real only when the expression 
G'^^{cOn) is symmetric with respect to interchange n <r^ h, that is only when cofi = a>n. 

To see how the procedure works for a larger number of oscillators, let us calculate 
the force acting on oscillator ni. First, we write 



■£G':;,{aj)M(^)e-"'', (6.10) 



as before. Now, we find ipn(<^) from Eq. (6.5 I where the force comes from the nonlinear 
term: 



MOJ) - —^ y r dt e"^' Unit) gnn' <An' (0 fn' (Ol = 

Zgnn' r dco' doj" 7 , ,. 7 . „. 7* . / , // x 
-— iA„(w)iA„'(oj )(A„/(w +aj -cu) 
, CO - (On J 2n 2n 

n' 

Expressing i^„, i^„', i^*, in terms of some other oscillators M2> «u «2> we obtain 

n,ni n' 



(6.11) 



(6.12) 



(72, "2 
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Again, inside the Green's functions one can replace a»' — > 6ti,i, , <J' — > w,!,, a>'+a>"-a> — > 
6ti„2. When they are taken out of the triple frequency integral, the latter evaluates to 
\p-nyit) ^n^{t) (A,*2(0- As a result, 

(a,-W,„)lA„i - 2] '^f.mn.-m^li'l'nx^ni^ (6-13) 



n2,'Jl>'J2 



'(2) 



where ^n „ „ ^ is the expression corresponding to the diagram for A^ = 2 in Fig. 
according to the rules of the previous subsection. Finally, the factor 1/2 has to be addec 
when deducing V'^^^ form TCa, in order to avoid double counting upon differentiation 

dV^'-^/drnr 

Generalization to the A^-oscillator case is straightforward: calculation of the force by 



iterative solution of Eq. ( 6.4 1 corresponds to iterations of diagrammatic equations ( 6.3a i, 
(6.3b I. The factor l/N in V'-^'' takes care of double counting upon differentiation. 

To conclude this subsection, we note that the perturbation theory, developed here, 
is closely related to that of Refs. lITTl |22]| . One obvious difference is that the linear 
part of the problem was assumed to be diagonalized in Refs. llT7ll22l . while here the 
linear coupling is treated perturbatively. This difference is purely technical; the present 
calculation corresponds just to the perturbative evaluation of the overlap sums which 
dress the nonlinear coupling in the approach of Refs. |[T7ll22l . Second, in Refs. lfT7ll22l 
one starts from a single excited normal mode as zero approximation and then calculates 
the force acting on other oscillators. Here, the independently excited oscillators are taken 
as the zero approximation, and the force acting on some given oscillator n is calculated. 
Clearly, the two approaches are closely related, so the structure of the perturbation theory 
is very similar in the two cases. It is the subsequent use of the information obtained from 
the perturbation theory which is quite different. 

6.4. On the number of diagrams 

It turns out that the total number of diagrams of the order g^~^, generated according 



to the rules given in Sec. 6.2 is given by the following exact expression: 



^^ = Z^^ Yj (3mi+2)(3m2 + 3)...(3mi_i+/c), (6.14) 

k=\ ' 0^mi^...'^mt-i^N-\-k 

to be derived below. Even though explicit, this expression is not straightforward to 
evaluate. Still, it can be shown that D^r grows exponentially with N. Let us estimate it 
from below and from above: 

Y 2^^W3A^-3-2A:)! y 2^"' (3A^ - 3 - 2A:)! {N-2)\ 

^ 2^(3A^-2-3;t)! ^ ^ ^ ^ 2^ ON -2- 3k)\ {k - \)\{N -\-k)\' 

(6.15) 
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The lower bound is obtained by simply taking one term of the m-sum with all m, = 
N -l-k; this is the largest term, so the upper bound is obtained by multiplying it by the 
total number of terms in the sum, where we use the formultj^ 



E 



1 = 



l^mi^...^mt^« 



{n- i+ky. 

k\{n-l)\ ' 



(6.16) 



At A'^ » 1 we approximate the factorials by Stirling's formula, denote k/N = x, and pass 
to the integral: 



1 

J 



dx (3 - 3x)3 



^/S^ V (3 - 2x)5 



3-2j 



2^--^(3 - 2x) 
x^X3 - 3x)3-3^ 



N 



<Dn< 



1 

/ 



(6.17) 



-■^r^3 



(3 - 3x) 



dx 



47rN^ V 41 - x)(3 - 2x)5 



2^-^(3 - 2x) 



3-2x 



(l-x)(l-x)x2^(3-3x)3-3^ 



The maxima of the expressions in the square brackets on the left and on the right are 
4.365. . . and 8.722. . . , respectively. For the numerical evaluation of Dn it is convenient 



to cast Eq. ( 6. 14 1 into a recursive form: 



£>A 



N~l 2/v-i 



l/t+1 



Z 

k=i 

n 

-z 



2^ 



-Sk 



k,N-l-k, 



m=0 



1 + 



3m 



* k,mt 



S l,n - 1- 



(6.18a) 
(6.18b) 



Numerical evaluation of the leading exponential gives 

^ -.(6.75.^.)x(l + i 

Dn w^oo \ A^ 



-1.5... 



(6.19) 



To derive Eq. (6.14i let us view the diagrams as trees formed by the dotted lines 
(Fig. [t]). The simplest one is a single-branch tree: in the order g^~' there are 2^~^ such 
diagrams. If there is one lateral branch, the number of such diagrams is 



^f = 9 Z f^^'i - 1^ + 2]2''-i2'2-i 6,,^!,. 



N~l , 



(6.20) 



h,h=i 



'This formula can be obtained by noting that the sum, denoted by 5(.„, satisfies the recurrence relation 
Sk.n = Sk-i.ii + Sk,„-i- Then one can (i) check directly that it is satisfied by ()j + ^ - l)\/[k\(n- 1)!], or (ii) in- 
troduce the characteristic function, xi'P, S) = Stn e''^'^*'^"^^''^ S k,n, for which the recurrence relation reduces 
to a simple algebraic equation, so \haXx{<l>, 0) = (I - e"^ - e'^)"', from which ij.,, is easily reconstructed. 
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Figure 7: Examples of two- and three-branch diagrams for effective 18-oscilIator coupUng. 

because 2''"' is the number of ways to build the main branch, 2'^"^ is the number of 
ways to build the lateral branch, 3(Zi - 1) -i- 2 is the number of electron lines where the 
lateral branch can be attached to the main one, and 1/2 stands because of the arbitrary 
choice of which branch is main and which one is lateral. Summing over the number of 
branches, we obtain 

N-l ^ 



k=\ h,...,k^O 



•)h+-+lk 



X2 ' - '' Sl^+,„+lt,N-i--k, 

which is equivalent to Eq. ( 6.14| ). 



(6.21) 



It turns out that Eq. (6.22 1 is equivalent to the following recursive relation, derived 
in Ref. 122113 






(6.22) 



The author was not able to prove the equivalence between Eq. (6.22 1 and Eq. (6.14) 



algebraically, however, it can be checked numerically for quite large N. For example, 
both equations give £)5() = 30426054945480277365983787382745806500. 



'"Strictly speaking, the quantity, for which Eq. 1 6.22 1 was derived in Ref. 1221 . is not defined in exactly 
the same way as Dn, since the perturbation theory developed in Ref. 1221 is not strictly identical to the 
present one. Still, the two perturbation theories are closely related, as discussed in the end of Sec.|6.3[ and 



the equivalence between Eq. 1 6.22 1 and Eq. 1 6.14 1 is a manifestation of this relation. 
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6.5. On the statistics of diagrams 

Let us define a resonance r as a sequence n''^ ^ n^ ^ ■ ■ . ^ n''^ ,n''^ ^ n^ ^ . . . ^ n''^ , 
where none of the m's coincides with any of the «'s. The freedom of interchanging n <r^ h 
can be fixed by requiring 

a>„'j + . . . + a)n'-^^ <a)n'^ + ... + Un'-^^ . (6.23) 

The frequency mismatch of the resonance is thus defined to be positive: 

mr^OJn\ + ■■■ + OJn'^^_ - {OJn\ + ■ ■ ■ + OJn'^). (6.24) 

To associate the resonance with a single site n, one can fix «j - n. In the case when 
some of the n's or some of the h\ coincide, one can introduce the multiplicity of each 
site n, mJI, defined as the number of times it occurs in the sequence of rfs or minus the 
number of times it occurs in the sequence of «'''s. For sites n that do not participate in 
the resonance, mj^ = 0. Due to action conservation, 

2]m;-0, 2]K| = 2A^,. (6.25) 

n n 

Viewing the sequence |m,'^) as a vector, we define the scalar product and the norm in the 
usual way: 

oo 

{m\m'')= YjKmn^ \m'\^ = {nf ,nf). (6.26) 

n=-oo 

The oscillators are coupled by an eff'ective Hamiltonian 

V = 2Vr COS {(prf^ + . . . + ^„^^ - (p-n^ "... - (prf^) , 

Vr - "Kr ^In[ ■ ■ ■ h"-^^ /,!', • • • hf^^ , (6.27) 

7C = - V %^^'> 

The summation is performed over {N,A)^/ Yin l"^nl' distinct permutations of ?ij, . . . , n'j^ 
and of «j, . . . , «]y . Besides the summation over permutations, TC- represents a sum over 



different diagrams whose number was estimated in Sec. 6.4 and for each diagram the 
summation over the positions of A^^ - 1 nonlinear vertices should be performed. Thus, 
TC",- is a sum of many terms, whose number will be denoted by N. Each term contains a 
product of L,. + Nr -2 frequency denominators, where L,- ^ A^,- is the minimal order in 
T which couples all oscillators together. 

It is convenient to explicitly separate all trivial deterministic factors and introduce 
the dimensionless random quantity 

2Lr+N,-2ANr-2 
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Figure 8: An example of the spatial structure for an TV = 9 diagram. 



Our primary interest will be the probability that 'Kr exceeds a certain value e 



-A. 



P\'Kr<e- 






(6.29) 



where 6{x) is the Heaviside step function, and the overline denotes the statistical averag- 
ing over the static frequency disorder. 

First, let us argue that the number of terms N is at most exponential in Nr, L^. It 
is convenient to introduce a pictorial representation of each term, which illustrates its 
spatial structure. Namely, let the direction from left to right represent the position on the 
chain. Then, each term of the perturbation theory, contributing to 7C, corresponds to def- 
inite locations of the free ends of the solid lines and definite locations of the dotted lines 
of the corresponding diagram. For the latter reason it is more convenient to represent the 
nonlinear vertex by a smal circle: 



(6.30) 



>Cr-X 




An example of a diagram with a spatial structure is shown in Fig. [8] The vertical di- 
rection on the picture has no particular meaning and is used just for the convenience of 
presentation. What really counts is the length of the projection of each solid line on the 
horizontal axis, as the sum of these projections for all lines determines the order in t. 

Now it is quite easy to see that out of a factorially large number of permutations of 
?i'''s and «'''s, very few actually contribute. Indeed, while permutation «i «-^ n2 in Fig. [8] 
does not change the order in t, the permutation n2 «-^ n^ already increases the order 
in T, and fii <r^ ng almost doubles the order. In fact, a typical permutation contributing 
to the factorial number, has a smallness t~^' ^ which suppresses even the large factorial 
factor. In fact, the number of permutations which preserve the order in t cannot exceed 
the number of vertices, i. e., Nr- I. 
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Figure 9: The structure of the diagram for a sparse resonance. 

Another source of different terms is the summation over the positions of the nonlin- 
ear vertices. In Fig. l8j for example, n'^ can be varied between fii and n2 without chang- 



ing the order in t. Clearly, this freedom is most important when the sites n^, . . . nj^ and 
are sparse, i. e., L,- » Nr (which will be called a sparse resonance). In this case 



the most important diagrams have the structure shown in Fig. [9} the solid lines should 
be drawn as vertically as possible, but at least one should go across the whole spatial 
region in order to connect all the oscillators. Moreover, for sparse resonances having 
|m„| > 1 clearly represents a disadvantage since imposing constraints on the oscillator 
positions decreases the freedom in making TUr small, without any gain in t. Then, each 
nonlinear vertex can be moved within a segment ~ Lr/Nr, so all vertices together give a 
factor ~ {LrlNr)^'. Using the fact that lim(D/v+i/DA,) ^ 6.75 . . . < e^, we arrive at the 
following bounds: 



li^N < 



Nr 



Nr 



,2N,- 



(6.31) 



Let us now turn to the statistics of of each individual term contributing to 'Kr- The de- 
nominators of the 2Nr external Green's functions contain differences of two frequencies 
only, whose probability distribution is 



m - |a»„ 



(Dn 



Pi -l(ro-) = - 



1-^ 

A, 



(6.32) 



where we denoted [xj = x6{x), i. e., [x] = ;c for jc > and [x] = otherwise. The 
Nr - 2 internal Green's functions contain linear combinations of many frequencies, so 
their probability distribution can be estimated from the central limit theorem: 



m ■ 



Z 



ninOJ,, 



Pm{TU) ^ - 



n\my 



exp - 



6^2 



\m 



|2A2;' 



(6.33) 
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since the dispersion of a single frequency is 

A/2 



doj T A" 



.2 



J A ' 12 

-A/2 



o;^ = — . (6.34) 



When the number of terms N is very large, the distribution at cr — > is the most impor- 
tant (the main contribution to "X",- coming from largest terms, thus small denominators), 
so we replace pim) by a box distribution of the height p{m -^ 0). The factor (2/A) has 
been included explicitly in Eq. ( |6.28 ), and for external Green's functions nothing else 



is needed. For internal ones, however, ^6/(;r|nip) should be included. Since the con- 
tribution of the internal Green's functions to L^ can vary from to 100%, and internal 
Green's functions can involve up to 2Nr frequencies, so that |mp can approach 2A^^, we 
can represent each term contributing to "XV as 

'X-^^'^'^^-a^, , (6.35) 

where the non-random coefficient a is bounded by 

< a < 1, (6.36) 



and ^ is a random variable constructed as 

^ ^- , (6.37) 

Xl . .. X^.^M^-2 

xi, . . . , xi^+N,.-2 being random variables uniformly distributed between and 1. 

Strictly speaking, ;c,'s are not statistically independent, since the same oscillator fre- 
quency can enter several factors in the denominator. The resulting correlations can lead 
to both bunching (meaning that if one of x,'s is small, than some others are also likely to 
be small) and antibunching (if one of x,'s is small, than some others are less likely to be 
small) of x,'s. First, we will estimate "X",. assuming all x,'s to be statistically independent. 
After that, we will study the effect of both type of correlations on the obtained estimate. 

The probability distribution for ^, defined by Eq ( |6.37 ) with independent x;'s. 






is discussed in detail in Appendix C It is crucial that all moments of this distribution 
diverge, so for the sum of N such random variables with random signs, ±^i ± . . . ± ^yv. 
the central limit theorem does not hold. Unable to calculate explicitly the probability 
distribution of this sum, we estimate the sum from above and from below. This corre- 
sponds to estimating the probability for the absolute value of the sum to be smaller than 
some value z, P {\ ± ^i ± ■ ■ ■ ± ^n\ < z), from below and from above, respectively. 
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As discussed in Appendix C the lower bound for the probability is obtained if one 
replaces the sum by its largest term: 



^{| ± ^1 ± . . . ± ^^1 < z) ^ n^,, . . . ,^^ < z) « exp 



-N I pn,-+l,- 



iiOd( 



(6.39) 



When yV is sufficiently large (exponential in A^,-, L^), it can be approximated as 

1 A^ In^'-^^'-^z. 



P{^u---,^N<z.}^e^> 



z (M- + Z., -3)!ln(z/e^'-+^'-3) 



(6.40) 



This probability is (i) very strongly suppressed for z smaller than a typical value at which 
the argument of the exponential is of the order of unity, and (ii) slowly approaches unity 
for larger z, so that all moments of z diverge. This typical value of z is found from 



1 



In^^+^'-^z 



z(A^,- + L,-3)!ln(z/e^'+^'-3) 
At large Nr, Lr the solution of this equation is given by 



1 

77' 



z 



[efiN'l^^r^'^-^ 



(6.41) 



(6.42) 



where we omitted the prefactor, a power of A^,- + L,, and the function f{x) is defined as 
the solution of the equation 

In particular, if W ~ c^'"^^' , then the typical z ~ {ef{c))^'''^^' . To give an idea of the 
behavior of /(;c), we note that /(I) - 1, and xln(ex) ^ f{x) < 2xln(e.j;) for ;c > 1. 
To obtain the approximate expression in the region around the typical value of z, we 
linearize the logarithm of the argument of the exponential in Eq. ( 6.40[ ) with respect to 
In z. Using the right inequality ( |6.31 1 the inequality 

^NI{L+N)J.NI{L+N)\ ^ „2 



and the right inequality ( |6.36[ ) we arrive at the following bound: 

^[Kr < e-^] > exp [-e-'+^CL.+iv,)] 



(6.44) 



(6.45) 



A lower bound for | ± ^i ± . . . ± ^yvl is obtained by "squeezing" its probability distribu- 



tion towards smaller values. As discussed in Appendix C this is done by first replacing 



distribution (6.38 1 for each term ^ by the log-normal distribution, applying the central 
limit theorem to the latter, and squeezing the resulting Gaussian into a box of the same 
height. This corresponds to assuming the sum to be uniformly distributed in the range 



< I ± ^1 ± 
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(6.46) 



where the left inequality (6.31 1 was used. Using the left inequality (6.361, we translate 
it into the uniform distribution of "X",- in the range 



o<n<r< 



n e 



2{L,.+Nr-3) 



(6.47) 



2 (y/JT^Nr)^^ 
Squeezing this range further, down to l/N^'', for pure analytical convenience, we obtain 

P ['kr <e-'^]< e^' '" ^'"-^ e{A - Lr In N,-) + OiU In N,- -A). (6.48) 

Having found the bounds for the statistics of "XV, let us obtain bounds for the follow- 
ing average: 

A 

\1 + ln"X",J^ - - f(A- Af ^r ['kr < e-^] dA, UJ = x e{x). (6.49) 



For the bound (6.45 I the integral can be calculated as 

A. oo 

- fa - A)"^ ^ exp [-e'-'') dA^ f i 



''e'-''-'dA 



max{0,i-/lr) (6.50) 

{A-Arf^^lCA-Ar-N), A-Ar-N^ VA^, 
A^! e^--^'-, N-{A- Ar) » V^V, 

where the abrupt cutoff of exp(-e''"'*'^) at A> Ar was simply replaced by a step function. 



For the bound ( 6.48 1 we obtain the integral 



\N+l 



{A-Ary'^-'liA-Ar-N), A^»l, 



OiA - Ar) I iA - ^)V^- dA^r^_ XJ. -^Vi- Ar. ^'-''^ 

In the following, we will be interested in the situation A - Ar » A'^, so our working 
estimate will be 



\_A - Lr In Nrf < L^ + In Krf <\A + 3{Lr + Nr)f. (6.52) 

Let us now see how the above estimates are modified by correlations between dif- 
ferent A','s in Eq. (6.37 1. To estimate the effect of bunching, we consider its extreme 
case by requiring that variables in groups of m are identical: xi = X2 = ■ . . = Xm, 
jCm+i = . . . = X2m, stc, This Corresponds to replacement ^ —^ ^i/m^ Lr + Nr - 2 —> 
(Lr + Nr - 2)/m is Eq. ( 6.37| ). Repeating the same steps, we obtain Eq. (6.45 1 with 



„A+3(L,+N,) 



gA/m+3(_Lr+N,-)/m 



, and Eq. ( |6.47| ) with e2(L.+A',-3) ^ g2m(L,+iv,-3) ^^^^ ^^^^ 
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m cannot be too large, as this imposes strong restrictions on the structure of the reso- 



nance, so we assume m <^ Lr,Nr- Thus, the modification of Eq. (6.45 1 has no effect 



on the estmate (6.52) since one still can replace exp(-e^ '^'''") by a step function in 
Eq. (6.50 1. As for modification of Eq. (6.47i, it is completely absorbed in the passage 
from Eq. ( 6.47| ) to Eq. (6.48 1. To estimate the effect of antibunching, we set some of 
the variables in Eq. (6.37 1 to x,- - 1, thereby replacing L^ + A^r - 2 by a smaller num- 
ber N' . This replacement can only decrease the value of ^, so the probabilistic upper 
bound on ^, expressed by Eq. (6.45 1, remains valid. In Eq. ( 6.47| ) we have to r eplac e 



,2(Lr+Nr-^) 



,2(iV'-l) 



, so this modification is also absorbed in the passage to Eq. (6.48 1. 



To summarize, the bounds established previously for uncorrelated x,'s, are loose enough, 
so that the effects of correlations fit in these bounds as well. 



6.6. On the number of resonances 

Here we estimate the number of different coupling terms of the form ( 6. 1 1, HiNr, Lr), 
corresponding to given values of A^^ and L,-. One of the oscillators is assumed to be fixed; 
otherwise, HiNr, Lr) would be infinite due to the translational degree of freedom. As will 
be seen shortly, HiNr, Lr) at large A^^, Lr can be represented as 



^(A^,, Lr) = 'RoiNr, Lr) [R{Lr/Nr)] 



Nr 



(6.53) 



where "RoiNr, Lr) is a power-law prefactor. Thus, most interesting for us is the quantity 



R{p)= Urn [<R{Nr,pNr)] 



l/Nr 



(6.54) 



Finiteness of this quantity (i. e., the fact that it is neither zero, nor infinity) is equivalent 
to the validity of representation ( 6.53[ ). The quantity p = Lr/Nr is nothing but the ratio 
between the order of the perturbation theory in t and the order in g. Note that L,- ^ Nr, 
since «^ t n[, for any k, k', and the oscillators must be coupled. Thus, 



R{p < 1) - 0. 



(6.55) 



To obtain an upper bound on R{p), let us choose some 2A^r-oscillator diagram and 
fix one of its free ends to a given site. Let us now count all possible positions of all other 
free ends and nonlinear vertices which produce the required L,-. In total, the diagram has 
2A^,- - 2 Green's functions, each corresponding to a certain displacement along the chain. 
Thus, we simply have to sum over all possible displacements, keeping their sum equal 
to Lr. As mentioned above, A^,. displacements cannot be smaller than unity, others can 
be zero: 



/ j / j ^kl+...+k3Nr-2,Lr 

kl,...,kN,^l kN, + l....,k3N,-2^0 



/ i ^ki+...+k3N^-2,Lr-Nr 

ki,...,k3Ni.-2^0 



2n 



-S 



dip 



-i(L,-N,:)<i, 



In (1 - e'0)3W,--2 



(Lr-rlNr-y)\ 
iLr-Nr)l{3Nr-3)\' 
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(6.56) 



This expression should be multiplied by 2^^^""^, since each kj can correspond to a step to 
the left or to the right. This procedure produces all resonances with the given L,. which 
involve the starting site. Moreover, each resonance is counted several times, because the 
spatial arrangements of the diagram, differing by the positions of the nonlinear vertices, 
but with the same positions of the ends, are counted as different configurations. Thus, it 
indeed produces an upper bound, and we have 



8(0 + 2)P+2 

Rip) ^ -^ — r ^ 8(2p 

21{p - l)/'-i 



1)^ 



(6.57) 



To obtain a lower bound, consider a resonance with nt = 2k - I, n^ = 2k, k = 
1, . . . , Nr- Coupling the oscillators by a diagram similar to the one shown in Fig.[9]gives 
Lr = 3Nr - 2. If one exchanges Hk <-> n^+i for some k, L^ is not changed, but a different 
resonance is obtained. Thus, we can construct explicitly 2^'~^ resonances with L, = 
3A^,- -2, which means that /?(/? = 3) ^ 2. One can generalize this estimate by considering 
resonances with%, -ni = L-1, and fixed |«)t-«/t+il = 1- ThenL^ = L+Nr-2, while the 
number of such resonances is atleast 2^'^'(L/2- iy./[{Nr- l)\{L/2-Nr)\] (Lis divided 
by two, since we are counting the ways to distribute nearest-neighbor pairs). This gives 



R(p ^ 3) ^ 



ip - yy-' 

\p-3)P-^ 



^ e{p - 3). 



(6.58) 



Taking into account the power-law prefactors, we obtain the upper and lower bounds 



as 



oNr 



^l2nN'r 



Nr 



-3 



Nr-l 



^ HiNr, Lr) ^ 64 



Nr 



hi 

Nr 



3Nr 



(6.59) 



where L-^J = x9{x), i. e., [xj = x for x > and [xj = otherwise. We will also need 
an upper bound for the contribution to "RiNr, Lr) from resonances with typical values of 
multiplicities |mjj| ~ m^. If we constrain \m'^\ ^ m,-, the number of oscillators is m'" times 
smaller. Thus, in the number of combinations we replace Nr -^ Nr/nir (but not in the 
ratio Lr/Nr, since it cannot be too small even for high multiplicities): 



Nr 



'R(Nr,Lr)^ 2] 



tNr/ltlr 



Lr 1 

A^~ 2 



Wrlnir 



(6.60) 



6. 7. On resonant triples for not too weak disorder 

The analysis of the stochastic phase volume of nearest-neighbor triples was per- 
formed in Sec. [5] under the condition t <^ p. This condition guaranteed that for two 
neighboring oscillators the phase space had a separatrix when the difference between 
their bare frequencies |cl»„ - (x>„+i| was sufficiently small. In the opposite case, t :s> p, 
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the minimal splitting of frequencies of the neighboring oscillators is 2tA because of the 
coupling (also known as level repulsion, in the quantum-mechanical language), so the 
nonlinearity is too weak to produce the separatrix, as seen from Eq. (twooscboundary=). 

Two normal modes of the disordered Unear chain can have eigenfrequencies much 
closer than tA only if they are well separated in space. Indeed, the effective coupling 
between two oscillators separated by a distance I can be obtained in the (/ - l)st order 
of the perturbation theory, and is given by Eq. ( |6.9| ). Strictly speaking, it is a random 
quantity, so its statistics should be analyzed. However, we note that the typical value of 
the splitting scales as t'A, rapidly falling off with the distance. As soon as this splitting 
becomes ~ pA or smaller, the separatrix can appear. On the one hand, this can always be 
satisfied when / is large. On the other, if I is too large, the coupling will be very weak, 
so the effective pendulum frequency will be very small, and the stochastic layer will be 
very thin due to Melnikov- Arnold exponential. Thus, the optimal value of I is such that 
r' ~ p. The same arguments can be used for the choice the third oscillator, needed to 
destroy the separatrix and to produce the stochastic layer. Thus, we face the situation 
similar to that of Sec. [s] but the coupling constant r should be replaced by t ^ t' ~ p. 
Hence, the average density of the chaotic spots for t » p is w ~ p^. 

Finally, we note that the situation described above (effective tunnelling coupling 
between two oscillators being of the same order as the nonlinearity and their frequency 
mismatch) corresponds to nonlinear coupling between two double-humped eigenstates 
of the linear problem. The distance between the two humps is given by /, and the two 
eigenstates correspond to different mutual signs of the humps of the wavefunction. The 
importance of double-humped states for the wave-packet spreading was discussed in 
Ref. m\. 



7. Probability distribution of the chaotic fraction 

7. 1. General remarks on the procedure 

Each site n of the chain can be characterized by a variable w„, the chaotic fraction 
of the thermally-weighted phase volume, summed over all possible guiding resonances 



involving the site n and sites to the right of it [to be formally defined below, Eq. (7.8 1]. 
For a resonant nearest-neighbor triple at t ■« p this quantity was defined in Eq. (5.29) 
and estimated in Sec. 5.4 the result being [see Eqs. ( |5.42a i, (5.42b i] 



Wn e 

p 






, IgEn^l max (cok-^f -y^iaJk-nf- (7.1) 

k=n,n+[,n+2 '—' 

k=n 



Note that the definition of w„ involves the thermal averaging but not the disorder 
averaging. Thus, w„ is still a random quantity, determined by random frequencies of the 
oscillators in the vicinity of the site n. Our goal in this section is to find its probability 
distribution function for values w much smaller than its disorder-averaged value, w <^ 
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min{Tp,p^}. For this purpose Eq. (7.1 1 is of no use, as it greatly underestimates the 
typical value of w (the value where the probability distribution has a maximum), giving 
w ~ {T/p)e~^^P. In fact, one can typically find chaotic phase volume at lower energies by 
taking into consideration guiding and layer resonances which involve many oscillators. 
Typically, for a given guiding resonance g and a given layer resonance {, the corre- 
sponding contribution to the chaotic fraction, w^^'^\ is very small (much smaller than for 
a nearest-neighbor triple, as high orders of perturbation theory have to be involved). That 
is, the probability that it exceeds some value wq - e~'^, 'Pw^'^^ > wq\, is small unless 
Wo itself is very small. However, w < wq means that w'-^'^ < wq for all of the numerous 
candidates for guiding and layer resonances. The probability of this event can be written 
as (assuming statistical independence of the contributions from different resonances) 



^{w<wo)-]~[(l-!P{ 



w'^s,n > ^ijU 



exp 



g/ 



Z^l 



«/ 



W^-S'^^ > Wo 



(7.2) 



Thus, S is determined by the competition of the smallness of each individual probability 
<p Ixv^s^^) > ^Q j and the large number of terms in the sum. 



It is more convenient technically to perform the summation in Eq. (7.2 1 in two steps. 
First, we calculate 



,iA) = Y,eiw^ 



igf) 



'), 



(7.3a) 



where 0{x) is the Heaviside step function, and the overline denotes the statistical averag- 
ing over the frequencies involved in the layer resonance, while the guiding resonance is 



kept fixed. Thus, e s is the probability of w 



(g) 



< e 



for a given guiding resonance, but 



after the optimization over the layer resonances. On the second step, we calculate 



S{A) = Y,'P{ 



g 



w^s) > g-^ 



= Z 



-Ss(A) 



(7.3b) 



where the averaging is performed over the frequencies of the guiding resonance. This 
latter averaging is assumed to be independent from the former one. Note that expanding 
the exponential and simply averaging S g would produce a wrong result, since at large A 
(the ones we are interested in) the main contribution to the average would come from 
configurations for which the expansion 1 - e~^^ ^ S g is not valid. This point will be 
discussed in more detail in Sec. 17.5.1] 

The main body of the bulky calculations is presented in Secs. |7.3f|7.5[ and the result 
is summarized in Sec. 



7.2. Definitions 



coupling as in Eq. (6.27 1: 



Consider a guiding resonance g (defined as in Sec. 6.5 r = g) with the effective 

V = -IVgCOsL.^ +--- + <Pnl - -^nf - ■ ■ ■ - -^4 )■ (7-4) 
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This coupling determines the effective pendulum Hamiltonian Hg and the frequency £i.g, 

4 = 1 l?«'l'p' - 2 Vg cos I ag = ^2\mS\^g\Vg\, (7.5) 

where = cp^s + . . . + 0„s -(p^s- . . .-cpf^s is the slow phase (pendulum coordinate), and 

p is the pendulum momentum (see Appendix D for the discussion of the kinetic energy 
term). 

The resonance condition reads as 



<jJy,S + . . . + (jjy,g 



OJr, 



I g + . . . + Is - Ij-S - . . .- l^g 



= /" 



(7.6) 



The right-hand side is positive, as we have imposed condition (6.23 1. Then the resonant 
configuration with the lowest activation barrier is obtained by setting to zero all /,-,'s, and 
all /,,'s except one, l^s, where nl is the site having the largest multiplicity mf , and if there 
are several such sites, the one with the smallest a>„ should be chosen. For this site one 
should set /,,« = /™"/m^g. This lowest activation barrier is given by 



7mm 

m „ I 

n% 



' rmin y 
J_ 



(7.7) 



„te) 



The chaotic fraction w)f corresponding to a given guiding resonance g is defined as 



w 



is) _ 



m^|V,(5lV + ...+V 



(7.8) 



where the thermal average is defined in the usual way, 
(■■■>r 



/(. . .)e^"-f^''«yT ^^_^ dUcpjln 



je("-l^''or)/TYl'^_^dInd4>n/2n 



(7.9) 



The pref actor |m^pW.v in front of the 5-f unction in Eq. (7.8) represents the effective 
thickness of the resonant hyperplane defined by the 5-function, and is discussed in detail 
in Appendix D Ws is the volume of the stochastic layer in the phase space of the effective 



pendulum. It is determined by the strength and the frequency mismatch of the layer 
resonance, so it depends on the actions of all oscillators involved in the layer resonance. 
The definitions and notations for the layer resonance are again as in Sec. [63} r = f. Then 
Ws can be found similarly to the case of a resonant pair, Eq. ( |5.19[ ): 

Qg 

7r|m«|2^2|Vgi 



Ws = 



■ A^.(A)ln ^^^^^^^^^ , 



(7.10a) 
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where the Melnikov-Arnold integral is defined in Eq. (5.13 1, and its index and the argu- 
ment are given by 



m = 2 



\{m^,mf)\ 
\ms\^ ' 



w„f + 



A 






. =£. ,7.10b) 



Definition (7.8 1 and expression (V.lOal are meaningful only when A » 1. When this 
is not the case, the two resonances should be treated on equal footing, so the situation is 
analogous to the case of the resonant triple. In particular, one arrives at a Hamiltonian 
similar to Eq. (5.281 with the kinetic energy given by {g/2)\rrf\^p^ + g{nf ,m^)p\p2 + 
{gl2)\n/^p^. In Eq. (7.8 1 the second (5-function should be introduced, and instead of 
|m^P Ws the prefactor should be {\rrf\^\mf\^ - (nf, mf)^)Wss, where W^s is defined analo- 
gously to Eq. ( 5 .30a i. It will be seen below, however, that for those values of w which are 
important for the transport, A » 1 holds. Also, A » m (indeed, m can grow at most lin- 
early with Ng, Ne, while A is exponential in Ng), so that the asymptotic expression (5.141 
for the Melnikov-Arnold integral can be used. 

Expression (7. 10a i for Ws is already quite cumbersome for explicit analytical calcu- 
lations. Hence, we will work with its upper and lower bounds. 



ION, 



A (gm 



g\ A^ I {2Ne) 



1 



l! 



■ e-"^!^ < \m 



^^^-W.<20NlH'^ 
« g \ A2 



I in 



(2A)' 



(m+l)/;n- 



-A/2 



[r(m)] 



\ln 



(7.11) 



They can be obtained using the following inequalities, which are easy to check: 

1 . |(m^,mOl 



;cln 


221/. 
X 


— < 


ag 


<A 


^ 



M'^ 



\v,\< 



2Ng ^ |m«|2 < 2NI, 
A^ 



2Nl 



< 



\mS\^ 



<Nt, 



The upper bound for Im^pW., can be further simplified by noting that large values of m 
are extremely inefficient. Indeed, if the guiding and the layer resonances have just one 
site in common, m ~ l/Ng. To increase m to unity and beyond, about mNg oscillators 
of the layer resonance should be restricted to overlap with the guiding resonance. Ac- 
cording to the estimate ( 6.59 1, such a restriction would decrease the number of available 
layer resonances by a factor ~ e"'^s , which would increase the typical smallest m/; by 



the same factor. Thus, we can estimate A„ 

mN„\(m+l)/n 



Aoe'"^s. 



The minimum of the function 



{Aoe'^^i)'' 



exp (-Aoe'"^V2) 



for Ao » 1 is reached at small m — > 0. Using also the fact that A^l'^e^^l^ < e~^^^, 
obtain the following upper bound: 



we 



im^rw, 



. A 
< 25NI - 
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gWe 

A2 



l/n 



-A/3 



(7.12) 



73. Thermal averaging 



To integrate over the actions, it is convenient to resolve the 5-function in Eq. (7.8) 
with respect to /„s . As will be seen later, the typical value of action for oscillators with 
n^ ,h^ ,n^ ,h^ + nl, as well as the typical thermal value of the action, r/|yu|, are much 
smaller than the typical value of 7^?^'": 



^n, . . ^^ iff ■ 



(7.13) 



Thus, we simply set I^s = 7" 



7.3.1. Lower bound 

In order to obtain the lower bound for w we write 



nA 



KTUt 



V^ 



1 1 r\mi\ 



nvji 



m e 12 ^ ^ rl'nnl/4 



■4f^,g%,(TIV\f-''J^lA^f"' 



n ^^. (7-14) 



while |Vf| can be averaged over the thermal distribution independently. For the integral 



over the oscillators of the guiding resonance we use the inequality, derived in Appendix 
|E]forA» 1: 



CO 

J 





exp 



-^n 



1 






ri e ^"dxn > exp -(2Ng - m^g + 4) ( — 1 



^.4/(2iVj-m,_g+4)- 



(7.15) 



and simplify the argument of the exponential using 

' 7T 1 



N4/(2A'j-m^g+4) 



{2Ng-m.+4)\ 



<2N„ A"^^''^-'""i^'UA'"'^, (7.16) 



while for the oscillators of the layer resonance we use the inequality 



CO 

J l]x^':"^"e-^"dxn>e-^^, 



(7.17) 



following from the fact that T(x + I) > e ". As a result, the lower bound for w is 



w^^''^ > 



n<t{Tle\iAf'- 
Ai2N{)l 



_ ^~E,IT ^-2N,kT' ^ Ar 



VJf 



^g'KgiT/liulf^ 



(7.18) 
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7. 3. 2. Upper bound 

Let us first average the actions of tiie oscillators belonging to the layer resonance, 
assuming it to have only one oscillator in common with the guiding resonance (see the 



discussion in the end of Sec. 7.2 1. 



w^^'"^ < 



25Nl 



HA 
r(l/(27r) + 1) gT 



Egg'KtiT 



T A2 \|^| 



Nc 



\ln 



e~^^^/^ {e-^'\ 



n^^- 



(7.19) 



As in Sec. 7.2 this estimate can be simplified using the fact that large values of \mj are 



inefficient. Indeed, suppose all \mj are either equal to m or (1 «; m < Ni). Then the 
product of the gamma functions can be estimated as 



[r(m/2;r+l)] 



2Nr/m (jn_\Nf/'^ 

\2nel 



At the same time, this leads to an increase of A, A ~ Km=\e^'e~^'l'". Clearly, the maxi- 
mum of the function ;7^^f/(2^)g-^i« '« '' '" js reached at the smallest m = 1. Analogously, 
one can show that restriction of one of the sites of the layer resonance to wf again does 
more harm by increasing A by a factor of the order of 1, than help by introducing a factor 
Eg/T. Using r(l/(2;r) + 1) < 1, we can simply replace all gamma functions by unity. 

The average (e~^^^) is sensitive to the multiplicities mf,. Again, large multiplicities 
are not efficient in lowering Eg, however, this fact is not obvious at this stage of the 



calculation, and it will be proven later, in Sec. 7.5.3 For the moment, we keep the 
typical multiplicity m^ as an additional degree of freedom and optimize with respect to 
it afterwards. The thermal average is estimated using the inequality 
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exp 
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{xi... Xn)" 



dx[ . . . dxN ^ A'^ 
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--(aA)i/(^'^+i) 
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proven m 



Appendix E Let seta = m'^ /4, N = {2Ng - m^g)/m^. In fact. 



(7.20) 
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where we denoted mi = m , for compactness. Thus, we arrive at the following upper 

yip 
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bound: 
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7.4. Summation over layer resonances 

The previous subsection was dedicated to the estimation of the contribution to the 
chaotic fraction w*^^'^^ due to a given pair of guiding and layer resonances (g, €). Here we 
are going to perform the summation over € in Eq. (|7.2|): 



SM) = Yj'PW'''^>(^^'} (7-22) 

Since afterwards this quantity will be summed over the guiding resonances, we are in- 
upper and lower bounds for w^^''^^ were 



7.3 



terested in those X for which 5^ «; 1. In Sec. 

given, w^' < w^^''^^ < w^f ' . Then, S ^ can be estimated as 






(7.23) 



where the averaging is performed over the frequencies involved in the layer resonance. 

7.4.1. Lower bound 

The lower bound, w< , depends on the frequencies through two quantities: the fre- 
quency mismatch ot^ and the coupling TCf . For the former we simply use the central limit 



theorem, Eq. (6.33 1. Moreover, since we are interested in very small me, the Gaussian 
exponential can be omitted. Then we obtain 



Sg> 



g'KJTiW 



A2 



'\NJA Zj N 



i^N.r^ , 



Ne 



A-^-\n 



j2m_ 

(p/e)^'- 



-In 



yNf-\ 



A^f-z-TCf 



NJ4 



(7.24) 



where p = gT/(\iu\A), and where [xj = x9{x), i. e., [xj = x for x > and [xj = 
otherwise. 
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To average over TCf , we use the left inequality in Eq. ( 6.52 1: 



S,> 



g^jT/\fi\rs 



A2 (2A^„) 



aNs/4 Zj Ne 



A -In 
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4Nj 
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2Nt 
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(7.25) 



We sum over L(, Ne using the left inequality (6.59 1: 
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(7.26) 



The sum over Li is estimated from below by its largest term, 

A^^ A- Eg/T - ln( yfhdfcp/4e) - Ne ln[32N^^ / (eT^ p)] 



Lf - 3Nf = 



NJ4 + N[ 



H2Ne/T) 



(7.27a) 
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V2^p Eg 32N^ 

A-ln^^—^-^-Ntln^ 

4e 1 ex^p 



\n{2NtlT) 

Nf+Ng/4 



Ne 



{Ne + A^„/4)^'+^s/4 



(7.27b) 



Again, we estimate the sum over A^^ from below by its largest term. It is easy to show 
that 



maxl--a| =gW«)^(«), N,nax^ 

N \n I 1 + 1/F(a) 

where the function F{a) is defined as the solution of the equation 

a 
In - = 1 + F. 
F 

We note that F{a ^ 0) ^ aje + 0{a^), F{e^) = 1, F(a -^ oo) ~ In a, and 

- Ina ^ F{a) < ln(2a), a ^ 1. 



(7.28a) 



(7.28b) 



(7.28c) 
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As will be seen later, the values of A which are important for the transport correspond to 
Ng,N[ satisfying the inequalities 



1 1 



N„<s: Nf <s: -, 



T p 



(7.29) 



Thus, with logarithmic precision we approximate 
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32Af- 1 

In — ^ — w In -7—, In 



4e 



er^p 



3-' "■ ^/2m~ep 



In 



(7.30) 



so we arrive at the following bound 
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1 g'K,{T/\^i\r^ ( ] 1 
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A + ln(\/p)-Ef,/T Ng 
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A + lnil/p)-E g/T Ng 

ln(l/T3p) 

eln(l/TV)^ 
ln(l/T) 



-1/2 



(7.31) 



7.4.2. Upper bound 

The procedure is not very different from that for the lower bound. After using the 
central limit theorem for rui and averaging over ^e with the help of the right inequality 
in Eq. (6.52i, we obtain 



S, <24 



F \mf/4 
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(7.32) 



We sum over Lc, Nc using the right inequality (6.59 1. The sum over Lc is estimated from 
above by its largest term, which corresponds to 

Ne 3N( n 



2 A^„/2 + 3A^^ ln(2e-3/T)' 



X 



i-^-i„l 



P 



50N, 



\6E, 
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- :^ln 



2n e^p^T 



multiplied by the number of terms. 
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This gives 



S^ < 



2mS \ T 



Z 

Ni 



mil A 



247T 



g'K,{T/\fi\rs 

' A2 

3Nf+l 



{2m«) 



g\N,l2 . 



\n{2e-^ It) j {6N[ + Nnf^l^+^Ne 



X 



(7.33) 



.-££-i„l 



P 



SON, 



16E, 



1/27T 



e 
8^ 



N„/m^^ 



Nt. 



2n e^Tp 



9-rn2 



Ngl2-¥?,Nt+l 



The sum over N[ is again estimated from above by its largest term multiplied by the 



number of terms using Eqs. ( 7.28a i-( 7.28c) and Eq. (7.291. Also, we estimate EgjT < 
Ngl{mlp) in the prefactor and under the logarithm, where we replaced the Gaussian 



distribution ( |6.33[ ) for the frequency difference in Eq. ( |7.6| ) by a box distribution. This 
gives 
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7.5. Summation over guiding resonances 
7.5.1. Toy calculation 

Before proceeding with the main calculation, let us illustrate the difference between 
two ways of calculation of S from 5 g, mentioned in Sec.|7.1[ namely. 



S = 



E: 



(7.35) 



and 



S - J](l-e-5.), 



(7.36) 



using a simplified "toy" expression, which nevertheless captures the essential features of 
the results obtained in Sec. 17.41 



(7.37) 



Sg^irPpf^.J^^^. 
' ^' ^\\n{\lTPp) 
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pln^- 



pX-pA^^ln'^ 



Figure 10; A schematic plot of e *s as a function of Xg, as determined by Eq. 1 7.371, 



Here Xg is a random variable uniformly distributed on ^ JCg ^ 1, ascribed to each 
guiding resonance. The number of resonances corresponding to a given Ng will be taken 
to be simply 2^« just for this toy calculation. 
The first calculation is straightforward: 



A'„=l 



J]s-g=f^ (IrPpf^ 1^(1 - Ap) [(e'^/'"(i/^^P) - l)pln ^ + 1 - ^P 
+ ^(^p- l)e'*/i"(i/-'W(l -e-i/pi"(i/-''P))pln-^| « 



+ 



(7.38) 



erPp^ In | e 

^ rPp 



AiMiin'p) 



The sum is dominated by A'^g = 1, which corresponds to resonant pairs. This means that 
only two resonances (i. e., 2^s as mentioned in the previous paragraph, taken a.tNg = 1) 
effectively contribute. Thus the whole calculation is meaningless, the correct one being 



-{\9\ 



w^g) < e--^ 



'W' 



(7.39) 
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It can be 



When 1 - e~^s «; 1, this is equivalent to Eq. (7.36 1. 

The plot of e~^^ as a function of Xg is schematically shown in Fig. 
viewed as a step function with a width of the order of pln(l/T''p). As long as this width 
is smaller than the offset, Xg = Ap - pNg ln^(l/T^p), defined as the point where S g - 1, 
we can simply approximate it by the Heaviside function. Then 



1-p 



A-Ng In^ 



1 

TPp 



(7.40) 
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(we remind that [x\ = xior x > and [xj = otherwise), which gives 






^A\n2/\nHl/T"p) 



(7.41) 



Note that the typical value of the offset, Ap - pNgln^{\/TPp) ~ pln^(l/T''p), which is 
indeed greater than the width, pln(l/T''p). 



7.5.2. Lower bound 



As discussed above, from estimate (7.31 1 for 5;, we pass to the estimate for S 



S >2][l-exp(-e'«<^)]. 
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(7.42) 



To perform the averaging, we replace exp(-e'^'') — > 6(-Ag), due to the abrupt cutoff at 
Ag > 0. To transform the lower bound on Ag into an upper bound on Eg, we note that 
the solution of the equation ax - In ^/x = A is smaller than (2/a)A for A » I, a ~ 1. 



This makes the restriction on £"„ tighter, so using the central limit theorem, Eq. ( 6.33 1, to 



average over Eg which for the lower bound can be approximated as Eg = Ifilrug/g, and 
the left inequality ( 6.52 1 to average over 'Kg, we obtain 
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(7.43) 



The sum over Lg, Ng is treated in full analogy with Sec. 7.4. 1[ to give 
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7. 5. 3. Upper bound 

Analogously to the lower bound, we write the upper bound for 5 in the form 



S < J][l-exp(-e''s)], 
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To perform the averaging, we again replace expC-e''*) — > 6i-Ag) and transform the 
lower bound on Ag into an upper bound on Eg, noting that the solution of the equation 
a;c + 2 In X = A is greater than A/ {2a) for A » I, a ~ 1. This makes the restriction on 



Eg looser, so using the central limit theorem, Eq. (6.33 1, to average over Eg which can 



be approximated as Eg = |/i|c7g/(mfg), and the right inequality (6.52i to average over 
iCo, we obtain 
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To sum over the resonances we use a generalized version of Eq. (6.60 1 to include mf : 
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The sum over Lg is estimated by the largest term multiplied by the number of terms: 
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Now we can optimize with respect to m^. The main competition is between the exponent 
~ Ng/m^, which is decreased upon increasing m^, and the term Ng In m^ which increases 
the base. The sum over Ng has a finite range, I ^ Ng s^ N'g"''. Note that Ng"^ is 
almost the same for m^ = 1 and m^ - Ng, since the term Ng \nm^ gives only a small 
correction to the term -(Ng/2) ln(l /rp^), due to Eq. (7.29 1. This means that when we fix 
1 ^ m^ ^ N^"^, and sum over Ng, analogously to the preceding subsections, the result 
of the summation has a form a^'^^'"' , where A^^ » 1 corresponds to the maximum of the 
whole expression and a > I. This maximum is determined by the balance of different 
terms which enter a, and thus also only weakly depends on m*. Thus, already as m^ is 
increased from 1 to 2, a^s » a'^s/^_ Hence, the sum is strongly dominated by m^ = 1. 

The situation with mf is less obvious, as the sum is effectively contributed by a few 
ml ~ I. To see this, let is fix ml and sum over A'^g. With logarithmic precision we 
approximate 
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The sum over Ng is treated in full analogy with Sec. 7.4.2 to give 
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(7.50) 



where the sum over mf converges atmf ~ I and can be estimated as 1/28. 



7.6. Summary 

We have arrived at the following expression for S (A) = -In "Plw < e~ 
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C = C(lnp/ Inr) is a slow function of the ratio of two logarithms. It can be bounded 
by estimating the argument of the exponential in Eq. (7.44 1 from below, and that in 
Eq. (7.501 from above: 
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[l+ln(l+x)]^ ^ 

4 ' ^ ^' ^ C(x) 



^3[l+ln(l+x)]^ 



(7.51b) 



C\ = Ci(ln(l/T), ln(l/p)) is a power-law function of ln(l/T) and ln(l/p) whose explicit 
form does not enter the final result within our precision. Note that the power p is in- 
herited directly from the value of the lower cutoff of L^, min{Lr/Nr), appearing in the 
expression (6.591 for the number of resonances, while the exponential asymptotics of 
this expression contributes to C. 

From Sec. 7.5 it is seen that the sum over the guiding resonances is mostly con- 
tributed by 

A E„ 



N, 



ln\l /rPp) 



A, 



(7.52) 



while the values of Lg are those close to the lower cutoff. Similarly, from Sec. 7.4 
obtain 



we 



N{ 



ln(l/TPp)' 



Ar ~ — 



NJ2 






A^.ln 



1 



TPpln\l/TPp) 



the values of L( again close to the lower cutoff. Then, from Sec. 7.3 we obtain 

,2/N, 



A ~ A^,A^/^^ 






(7.53) 



(7.54) 



The inequality A-Li In N[ ^^ Ng, necessary to use the first regime of Eq. ( 6.50 1, becomes 
evident at this point. As a result, the long calculations of Sec. |7]can be summarized in 
very simple terms: the chaotic fraction can be written as the product of three factors 



w 



S\XA)(e~^c/^,)e~Es/T 



l^e~^cin,y 



(7.55) 



where the pref actors containing T,p in powers ~ 1 are neglected, and the most important 
configurations determining to fiw < e~'^} correspond to the logarithms of each of the 
three factors being of the same order, -A. 
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It may be interesting to compare the contributions to the stochastic layer width from 
the layer resonance and from the next best candidate. For the latter the orders in r, g are 
Ukely to be the same, Lr = Lc, Nr = Lf, while the frequency mismatch, tu,-, is greater 
than that for the layer resonance, me, by a factor ~ 1 . Thus, we set mr = 2me for an 
estimate, so that A^ is twice greater for the second resonance, and focus on the ratio 
between the values of the Melnikov-Arnold exponentials, 

p;|l~exp[-(2^-"'.- da] ~. -'"=<""".. <7^56) 

Thus, the contribution of the next best resonance to the stochastic layer width is para- 
metrically smaller. 



F inally , to justify inequalities (7.29 1, Ng ^: Ne ^: 1/t, 1/p, as well as assump- 
tion (7.13 1, T/\iu\ «; /„g «; /™", we have to specify A. The typical values of A are 



those for which SiA) ~ \,i.e. A~ ln(l/p)ln^(l/T''p), so inequalities (7.29 1, (7.13) are 



indeed satisfied. As we will see in Sec. [9j the transport is dominated by regions with 



anomalously small w, i. e. with larger A, so that Eq. (7.5 la i works. It ceases to be valid 
for A < ln^(l /t^p), so it is not suitable to find the average value of w, which is domi- 



nated by rare resonant configurations (triples). Moreover, for /I » In (1/t''p) the typical 
strength of the guiding resonance V^, which scales as (t''p)'^«, is smaller than p even 
at p «; T, so that the range of validity of the pendulum approximation coincides with 



that of Eq. ( 7.5 la I. The relation between t and p becomes important in determining the 
nature of the rare resonant configurations (whether they correspond to nearest-neighbor 
triples, or spatially separated ones). 

8. Stochastic pump and Arnold diffusion 

The standard description of Arnold diffusion in multidimensional dynamical sys- 
tems is based on the stochastic pump model [54ji55J. Namely, motion of the oscillators 
participating in the guiding resonance (the chaotic spot) is assumed to have a stochas- 
tic component with a continuous frequency spectrum. The continuous spectrum arises 
because the motion corresponds to successive passages of the separtrix of the effective 



pendulum at random instants of time and is analyzed in detail in Sec. 8.1 These oscil- 
lators of the chaotic spot are coupled to all other oscillators in the system (typically, in a 
high order of the perturbation theory). The combination of the two ingredients (coupling 
and continuous spectrum) results in energy and action exchange between the chaotic spot 
and the surrounding oscillators. This exchange occurs by small random amounts, which 



corresponds to diffusion. The corresponding diffusion equation is written in Sec. 8.2 
and its consequences are analyzed in subsequent subsections. 
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8. 1. Stochastic pump spectrum 



Let us assume the general form of the guiding resonance, as defined in Sec. 6.5 

J]m^„cbn = 0, (8.1) 



where ioni^n) - (^n +gln, and the integer numbers {/nf ) satisfy Eq. ( 6.25 1. If the dynamics 
of the guiding resonance can be described by an effective pendulum Hamiltonian, and 
the stochastic layer around the separatrix is thin, for each oscillator participating in the 
guiding resonance its phase 0„ can be written as 



m° 



0„(O ^ OJnt + -—^ 0(0, 



|m«|2 



where 



0(O = _^mf0„(O 



(8.2) 



(8.3) 



is the slow phase which represents the coordinate of the effective pendulum (see Ap- 



pendix D|for details). For rotation and oscillation of the pendulum, corresponding to the 



pendulum energy being above or below the separatrix energy, respectively, ^{t) can be 
represented in the following form (see also Fig. Ill: 



rotation. 



■^ v^ J 4 arctan e"**^' 'j^ - n, 

^ .Zj 1 [4 arctan e"«('"'j) -n]- [4 arctan e"s('"'j''2-0+i/2) _ ;,;•], oscillation. 



J=~tX 



(8.4) 
Hi. The instants of time tj, which 
are determined by the separatrix mapping [Eqs. (5.16ai, (5.16bi with t - 0/(An^)], 



Here Q.g is the effective pendulum frequency, see Eq. (7.5 



represent the random ingredient of the problem. Assuming the typical interval tj^\ -tj » 
1/Qg, we can represent the Fourier transform as a sum of independent contributions: 






-i(ml2)^{f)^icot ^j. 



1 



j=-oo 






(8.5) 



The Melnikov-Arnold integral J{m{^) is defined in Eq. ( 5.13| ), and here we use the stan- 
dard notation for its index m - 2m^/\nf\^, which is not necessarily an integer. If the 
times tj are assumed random and uncorrected, the power spectrum is given by 



'(OJ) 



oo 



-i{m/2)^(t)+i(m'/2)^(f 



\ 



M'-t') 



d(t - t') = 



{6t ^)t { 2J[,n{oJIQ.g)J[„,'{ojlQ.g), 

Ji,„(co/Q.g):7l„f{aj/Q.g) + ^nt{-co/^g)^nf{-co/ng). 
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(8.6) 
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Figure 1 1 : Schematic representation of the behavior of the slow phase 4>i.t), corresponding to rotation (a) 
and oscillation (b) of the effective pendulum. The flat regions around odd multiples of n correspond to 
long intervals of time spent in the vicinity of the unstable equilibrium of the pendulum, while the steps 
correspond to passages through the stable equilibrium at random instants of time tj. 

Here {6r^)t is the average repetition rate of t/s for rotation, and the double repetition 
rate for oscillation. Note that {6r^)t ~ Qg, up to a logarithmic factor. The angular 
brackets here correspond to the average over the dynamics: 



{f{ti,...,tN)}t 






f{ti+t',...,tN + t'). 



(8.7) 



The assumption of random and uncorrected tj is not fully correct. Small residual 
correlations in the phase found in Ref. |[67l for the standard mapping translate in cor- 
relations between different t/s. These residual correlations are usually explained by 
the system sticking to various structures inside the chaotic region, such as stability is- 
lands where the diffusion is anomalously slow, or accelerator modes, where the diffusion 
is anomalously fast. When dealing with the standard mapping, one often introduces a 
small amount of external noise which prevents the system from sticking IISTIIFT'I . The 
present problem, however, is more complex than the standard mapping, as we already 
discussed in the end of Sec. |5.2[ In addition to that discussion, we note that the resonant 
triple is not isolated, but it is subject to perturbations from surrounding oscillators. Even 
though these perturbations are not truly random, it is quite plausible to assume that they 
are sufficient to prevent the triple from sticking, and thus effectively fulfill the role of the 
external noise. 
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As we have seen in Sec. |5.4[ a chaotic spot is not always described by an effective 
pendulum with a thin stochastic layer around the separatrix, and the guiding and the 
layer resonances are not always distinguishable, so that above expressions are not always 
applicable. Still, for any oscillator n, belonging to the chaotic spot (n e chaotic), we can 
write 

Mt) - VA^e-'""'^„(0, (8.8) 

where ^(t) is a certain random process. Its statistical properties cannot be determined 
precisely. Still, one can assume {^n{t))t = 0. Also, for any set of integer numbers 
m = {m„} the two-time correlator and the corresponding spectral density can be defined: 






nechaotic 



iai{ti -t2) 

(8.9) 



When the guiding and the layer resonances are not distinguishable, the chaotic dynamics 
has only one characteristic frequency scale, Q.g. Then, the leading exponential in the 
high-frequency asymptotics \oj\ » Q.g, can be estimated as 

vm(a;) ~ T^ e"'^''"'^^^"^^ (8.10) 



For a resonant triple, considered in Sec. 5.4 Q.g ~ ^/fpA, which is vaUd at r «; p As 



discussed in Sec. |6.7[ in the opposite limiting case, t <s: p, one should effectively replace 
T by p, so £lg ~ pA. 

8.2. Dijfusion equation 

Let us now consider the dynamics of oscillators, surrounding the chaotic spot. The 
latter is assumed to be described by the Hamiltonian 

H = Y, [^nln + I /') - 2VgCos 2] mM. (8.11) 

Other oscillators are coupled to those of the guiding resonance by the Hamiltonian of 
the form 



y--2]2Kcos 2]m>„ 



(8.12) 



where the integer vectors nf and m'' have a nonzero overlap. They satisfy the constraints 



imposed by the total action conservation, Eq. (6.251. Generally speaking, Vg and V,- 
depend on the actions of the participating oscillators. For the moment we assume the 
effective pendulum description of the guiding resonance to be valid, so we take Vg and 
Vr exactly at resonance (the opposite case is considered in the end of this subsection). 
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Let us restrict our attention to a finite number N > Ng oi oscillators. To separate the 
slow pliase of tiie guiding resonance, we perform a canonical transformation (see also 



Appendix Di: 



N 



N 



4>^Yj '""^" ~ *--^' ^-'' ^'k^^ "'^"^'^ ~ ^-k' t^' k^2,..., 



N, 



n=l 



n=\ 



(m^, a>) 



(8.13) 



In = m„P\+ 2jPkUkn, P ^ Pl+ I g|2 ' O) ^ (a;i,...,w/v). 



k=2 



g\mg\^ 



Here p is the momentum, conjugate to the slow phase, and counted from the point of the 
exact resonance. Ukn is some constant matrix satisfying the appropriate orthogonality 
relations: 



iUk>vlLk'>\> = ^kk', {Uk>vWi^) = 0' (" ^). 



nk 



mil\mS\^, k=\. 



(8.14) 



For oscillators not participating in the guiding resonance one can set Ukn = Stn, i- £-, 
Pk = h- It is also convenient to introduce the vector iJ with all components m^ = 1, 
which corresponds to the total action: 



N 



N 



hot = Yj^'' = Tj P'^^^k^ ^'^' ^^'' ^'^ = ^^''' ^'^ = ^- ^^-^^^ 



n=\ k=2 

In the absence of the couplings Vr the hamiltonian is given by 



H = Y_\{^,Uk)Pk + ^pl 

k=2^ ^ 



-, ~9 +Hi,, H„^^^^p^-2V„cos(f>, (8.16) 



so the actions p2, . . . , /?w, as well as the pendulum Hamiltonian Hg, are conserved. When 
the couplings Vr are switched on, the change in the actions between two moments f i < t2 
is determined from the Hamilton's equations, dpt/dt = -dV/dcp'^., and the change in Hg 
is found from dHg/dt = -{dV/d4>){dH/dpy. 



n 
Spk>i = -i2_jVr{m'',Ui^) I <i?'exp 



\mS\'^ 



+ c. c, (8.17a) 



dH, 



-'■Z 



Vr 



{nf, m^) 






exp 



(m'',mg) ^ , 



+ c. c, 
(8.17b) 
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where "c. c." stands for the complex conjugate and lUr = Z„ m'^(b„ = {nf, u). Note the 
relations 



N 



2_^ dpkiUf,, M ) = 0, 

k=2 

N 

J^6pk{ui,,ay) + 6H = V{0)-V{t). 



(8.18a) 
(8.18b) 



k=2 



The first relation follows from the orthogonality relations and represents the total action 
conservation. The second relation is obtained by integrating by parts and represents the 
total energy conservation. In fact, these relations hold separately for each given term r 
of the perturbation. 



The integrand in Eqs. (8. 17a I, (8.17bi is randomly oscillating due to the chaotic 
dynamics of ^(t), so (Spk) remains bounded even for large t2 - ti. In contrast, the 
average of the product {6pkSpt)t has a contribution growing linearly with time t2- ti: 



{5pk6pk')t ^ 2{t2 - ti)^V'jvrjj,-{mr){nf ,u^){m^ ,u^,), 



(8.19) 



where the stochastic pump spectral density is defined analogously to Eq. (8.6 1: 



v„Aoj)= I (e'Wi)-^te)]«,m^)/KP^^gMf,-f2)j(j^_j2). (8.20) 



The interference terms with r t r' give a contribution which does not grow with time 



and thus are neglected. Using Eq. ( 8. 1 8b i, we can also write 



(SpkSHg) 



--Z 



{SpkSpk'). 



h-t\ ^ ^2 - fi 



{u,,cb). 



(JSHgfh ^ y {6pk6pk')t 



{U,,UJ){U,,,0J). 



(8.21) 

To w rite the diffusion equation, it is convenient to pass from the pendulum energy 

Hg, Eq. (7.5 1, to the pendulum action W using the standard expression in terms of the 



complete elliptic integrals K{k) and E{k), 



W=^- ^^g f E{K) - (1 - k')K{k) k<\, ^ 



Hg + 2Vg 



nVg\mS\^ I iK/2)Eil/K) + 1/2 k>1, 



4K 



(8.22) 



Since transformation (8. 13 1 is canonical, the volume element in the phase space is pre- 
served: 



dWd^ Y-(dpkd(p[ 



N 

n 



In \\ In 

k=2 n=\ 
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In ' 



(8.23) 



Not being interested in the oscillating phases, we introduce the distribution function / 
in the action space, so that the probability to find the system in an infinitesimal volume 
element is given by 

dP = fiW,p2, . . .,Pn)dWdp2 . ..dpN. (8.24) 

It is the flatness of the phase space measure appearing here that was the main purpose of 
passing from Hg to W. 



Now, from Eqs. ( 8. 19 1 and ( 8.21 1 we can immediately pass to the diff'usion equation 



dt dW 



aW opk 



dpk 



r 

dW 



Dkw = Dwk = -Dkk'{uu,,u)) 



Diyiy = Dkk'iui^,aji){Uj^,,oj) 



dHg 

dW 
dH, 



. df . df 

^'^ dW ^ ^''' dpk] 


, (8.25a) 


|2(^r), 


(8.25b) 




(8.25c) 




(8.25d) 



where summation over repeating indices k,k' = 2, .. .,N is assumed. 

Diffusion in W is confined to the narrow stochastic layer of the width Ws, given by 



Eq. (7.10a). At times longer than the diffusion time across the layer we can assume the 



distribution function to be independent on W, and write the probability ( 8.24 1 as 

dP - f{p2, ..., Pn) W, dp2... dpN. (8.26) 



As discussed in detail in [Appendix F[ projection of the full diffusion equation on the 
layer leads to a closed equation for /: 



W, 



df 
dt 



Z 



d 



df 



7^^ OPk Opk 



(8.27) 



hk'=2 



which describes the Arnold diffusion in the space of the variables p2, . . - , pn- The diffu- 
sion coefficient is given by 



dW 



Dkk' = Dk'k = I 

layer 



{u,^,D,u,^,) 



(mj, D, aj){&, D, Uj^,) 
(o), D, d)) 



(8.28a) 
(8.28b) 



and we use obvious tensor notations, 

{u, [n/ (gi m^, m') = (m, nf){n^, u), {[n/ ® m'], m') = rr/{m' , u). 
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(8.29) 



As seen from Eq. ( 8.28a i, the diffusion coefficient vanishes if the tensor D is separable. 
This would be the case if the sum over r involved only one term. If the sum involves 
terms of very different strength, the contribution of the strongest one is eliminated. This 
corresponds to the situation discussed in Ref. 1541 : the strongest perturbation of the 
guiding resonance effectively determines the stochastic layer width Ws, but does not 
contribute to the Arnold diffusion. This strongest term is called the layer resonance, 
r = €. Also, the diffusion coefficient (|8.28a|) vanishes identically for the directions 



corresponding to the change of the total action and of the total energy: 



N 



N 



2_^ Dkk' (m^, ,j^) = 0, 2_j ^kk' (ut , ^) = 0. 



(8.30) 



k'=2 



k'=2 



Thus, the minimal number of oscillators, necessary to produce the Arnold diffusion, is 
A^ = 4 in a system with two conserved quantities. 



Constraints (8.301 ensure that any distribution function / which depends on pk only 
via the combinations 



N 



hotiPl, ...,Pn)^2_j P''^-k' ^^)' ^totiPl, ...,Pn)^2_jY-'' -k^P'' "•" 



k=2 



k=2 



■pI 



, (8.31) 



is a stationary solution of Eq. (8.27i. Two specific distributions are particularly impor- 
tant. One is the microcanonical distribution, determined by two parameters l'^g,,H^tor 



fi°,.H°SP^' ■■■^PN) = e ^6 [HtotiP2, ...,Pn)- ti^o) '5 [ltot{P2, ...,Pn)- Iwt) > 

(8.32a) 
where e~^ is the normalization factor, determined from 

/ - r <5 {hUP2, ...,Pn)- K) 6 (ltot{p2, ...,Pn)- t) W, dp2... dpN. (8.32b) 

Thus defined S = 5'(/°p,,//J)^,) is nothing else than the entropy. The microcanonical 
distribution is the one to which N oscillators relax after having started from an initial 
condition characterized by some values of Itot(P2, ■ ■ ■ , Pn) - I^jt ^^'^ HtotiP2, ■ ■ ■, Pn) - 
H^gj, provided that these N oscillators are isolated from the rest of the chain. Indeed, the 



microcanonical distribution ( 8.32a i is the only one which (i) is a stationary solution of 
the diffusion equation (|8.27|), and (ii) has definite values of Ifot and Htot- In other words. 



the diffusion equation makes any initial distribution evolve towards the flat one on the 



allowed manifold, which is precisely the microcanonical distribution (8.32a I. 

The other important distribution is the grand canonical distribution, determined by 
two parameters p = l/T (the inverse temperature) and fi (the chemical potential): 
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,Pn)] 



(8.33a) 



where e^^ is the normalization factor, determined from 

^-/3F ^ C^-ISlHUP2,■■■,PN)-^'I,o,iP2,■■■,PN)]^r^ dp2...dpN=\. (8.33b) 

Thus defined F = F(fi, /j.) is the free energy (more precisely, the grand canonical ther- 



modynamic potential). The grand canonical distribution ( 8.33a i is established on a finite 
segment of the chain, considered to be a part of a larger segment, upon relaxation of 
this larger segment to equilibrium. This larger segment may be considered isolated, 
then the grand canonical distribution for the smaller segment can be derived from the 
microcanonical distribution for the larger segment using the standard argumentation of 
statistical physics Il82ll . 

In the region of the phase space where the guiding and the layer resonances cannot 
be easily distinguished, the above considerations should be slightly modified. Let us 
consider a resonant triple occupying the sites n = 1,2,3, and other oscillators on the 
sites n = 4,. . .,N. As discussed in Sec. 5.4 the internal chaotic dynamics of the triple 
cannot be described in the diffusion approximation. The relevant information about this 
dynamics is contained in the spectral density v,„r{aj), defined in Eq. <\S.9\. Besides, 



the isolated triple can be characterized by the values of its total action Itr and the total 
energy, which are conserved. The Hamiltonian of the triple is given by Eq. ( |5.28 1, which 
we write as //tr('^tr) + H. 

Under the action of the perturbation, /„, /tr, H acquire the increments 



6In>3 = -i^j ^X fdt'e''^'''E^.{t' 



) -i-c. c. 



(8.34a) 



J] Sin + 6hr = 0, 



n=4 

N 



2 OJnSIn + ^tr^/tr + SH = V(0) - Vit), 



11=4 



(8.34b) 



(8.34c) 



where cotr = dHtddhx, and Emi-it) was defined in Eq. (8.9 1. As discussed in Sec. 5.4 the 



excursion in H is limited to the chaotic region of a finite phase volume, so the Arnold 
diffusion occurs in the space of (/tr, I4, I5, ■ ■ ■ , In)- Thus, we introduce the distribution 
function /, such that the probability for the system to be in an elementary volume of this 
space is 

dP - /(/tr, /4, . . . , In) W,, dh, dh... diN, (8.35) 

where Wss is the chaotic phase volume in the space (/i, /i, h) at fixed /tr = I\ + I2 + h, 
as defined in Eq. (|5.30a[). If we let the index k assume the values tr, 4, 5, ... A'^ and define 



the vectors m^ as 



Ufcn 



5\n + S2n + S^n, k = tr, 

Skn, k = 4,...,N, 
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(8.36) 



then the diffusion equation for distribution equation has the same form as Eq. (8.27) 
with the only replacement W^ -^ W^s, and the diffusion coefficient is again given by 



Eq. (8.28a) 



8.3. Relaxation of remote oscillators 

Here we apply the diffusion equation from Sec. |8.2[ Eq. (8.27), to study thermaliza- 
tion of an oscillator which is sufficiently far away from the chaotic spot. In particular, 
we find the thermalization time. 

Let the Mh oscillator have the weakest coupling to the chaotic spot, i. e., D^n, Dj^i^ «; 
Dkk' for k, k' < N. Then the first N - \ oscillators equilibrate among themselves much 
faster than with the Mh one. In this situation, equilibration of the Mh oscillator can be 
described by seeking the distribution function in the form 



f{p, In; t) = f^M f\h; t) + 6f{p, In; t), 



(8.37) 



where we have denoted p = {p2, . . . ,pN-\)- Being sufficiently far away, the Mh oscil- 
lator is not involved in the guiding or layer resonances, so pN = In, and Ws is assumed 
to be independent of In. In the first term, f^ is the grand canonical distribution func- 



tion, Eq. (8.33a), for the first A'^ - 1 oscillators. Finally, 6f is the small correction, 
0{DNN/Dkk'), whose component along f,f in the space of functions of p is fixed by 



/ 



6f Ws dp2 . . ■ dpN-i = 0. 



(8.38) 



Neglecting terms of the second order in Dnn and using constraints (8.30 1, we write 



Eq. ^JJ} as 



WJs 



P Qt 



-I- 



k,k'=2 



d_ 
d 



+ /3{cbN - fJ^) 
d6f 



f+fL'^sD. 



I}.m' 



'NN 



dK 



+ /3{cbN -fi) 



„ W,Dkk' „ , 
opk opt 



where we denoted 



9iV 



N-\ 



k=i 



d{W,DkN) 
dpk 



/"+ 



(8.39) 



(8.40) 



Note that the last term in Eq. ( 8.39 1 is of the same order as others; its role is to compen- 
sate the dependence of Ws, bN, Dnn on p. However, being a total derivative, it vanishes 
when integrated over p2, ■ ■ . ,Pn-i- Thus, integrating Eq. (8.39l over p2, . . . ,Pn-\, de- 
noting 

<...>x= f{...)f^^Wsdp2...dpN-U (8.41) 
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and noting that 



J 



On ffi^i, dp2... dpN-i = — ^, /3ia)N - mDNN)±, 



dh 



(8.42) 



we obtain a closed equation for /"(In)'- 

Of d 



dt dl 



(Dnn)± 



N 



dh 



+ Pi(bN-fi) 



f 



(8.43) 



Solution of this equation is not trivial since both {Di\in)j_ and a>i\i depend on If^: con = 
CON + g^N, and {Dmn)± '^ In, as discussed in Sec. 8.4 Still, its stationary solution is 
readily found to be given by ^,;j(/a?). 

An analytical solution of Eq. ( 8.43 1 describing the relaxation dynamics can be ob- 
tained when gl^ ■« con - y". so we replace ojn -^ ojn- After reseating the action and 
time variables, we can rewrite Eq. ( 8.43| ) as 



dfU = dx [x{dx + 1)m] , 

X = /3(a)!\i - iu)In, t' = /3{a>N - y) — —. t = Yt. 



(8.44) 



dh 



Seeking the solutions oc e , and performing the Laplace transform. 



I' 



u{s) = I u{x) e "'^ dx, 





we obtain 



5 -— [(5 + 1)m] - Ail 
ds 



ua{s) = 



{S+ 1)-^+!" 

The inverse Laplace transform takes especially simple form when A is integer: 

1 d-{x-e-^) , 

= -^_ — -ni — = ^nix) e , 



(8.45) 



(8.46) 



(8.47) 



where L„(;c) are the Laguerre polynomials, which form a complete and orthonormal set 
of functions on < jc < 00 with the weight e~^: 

/oo 
Ln(x)Ln'(x) e'"' dx = 6nn', 2^ L„(xo)L„(x) ^"^ = S(x - Xq), x,xo > 0. (8.48) 



The eigenfunction uq(x) = e ^ with zero relaxation rate corresponds to the equilibrium 
distribution, 

d (In) = /3iajN - y") e-^^'^""^'^". (8.49) 
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The eigenfunction ui(x) = (l-x)e ^ corresponds to a small mismatch in the temperature 
[or in the chemical potential, as they enter in one combination, f^ioJN - fJ-)] '■ 

g Sp = 6fi{caN-li)[l- PicoN - I^)In] e-P^'^'-i'^'' . (8.50) 

Since the components of the distribution function, proportional to eigenfunctions m„>o, 
decay as e~"^' , the typical action (or energy) relaxation time of the oscillator can be 
estimated as r^^ ~ {T l\n\f{DNN)'^ ■ 

From Eq. ( 8.43[ ) one can easily obtain the fluctuation-dissipation theorem - a relation 
between dynamical fluctuations of the system in equilibrium and the system's dynamical 
response to an external perturbing force 1831 . First, let us find the power spectrum of the 
fluctuations of the action, 

oo 

SN{i^) = ^[{iN{t)iNm - {hfy'^ut. (8.51) 

— QO 

The average can be written as 

</;v(0 hm = f /^4/4) < r In G{In, f|4) diN, (8.52) 

where the first integral represents the averaging over the initial conditions 7^(0) = 7^, 
and in the second integral we introduced the conditional probability G(In, t\I^) - the 
probability distribution of In at time t, given t hat at t = it had a fixed value In = 7^. 
Thus, G{In, t\I^) is given by the solution of Eq. ( 8.43 1 with the initial condtition /"(7a?) = 
6(In - 7^). This solution is easily found using the second of the equations ( 8.48 1: 



G(In, f|7° ) = /3(0JN -fi)Yu L«CS(ww - fi)ll) U(/3{aJN - fi)lN) e-^^""^-^^'^-"^' . (8.53) 

n=0 

This immediately gives 

1 + e-^l'l 1 2 

Let us now find the susceptibility XNi(^) which determines the response of In to a 
periodic perturbation of the system's Hamiltonian, 

<577 = - {F^e-"'' + Fle"^') In, (8.55) 

where F^i is the generalized force corresponding to In- It has the dimensionality of 
a frequency, so it can be viewed as a periodic modulation of the frequency ojn- This 
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modulation produces a periodic deviation of the action In from its equilibrium average 
value; by definition, 



(In) - (In)\f^=o + XNioj) F^e-"'' +Xn(^)FJ''' + 0{FI). 



(8.56) 



The average of I^ on the left-hand side is performed over the perturbed distribution 
function, 



fVN) = fl/h) + 5fl{lN) e-'^' + {6fl{lN)T ^'"' + 0{Fi\ 



(8.57) 



which should be found from the perturbed Eq. ( |8.43[ ). The frequency oj^ enters Eq. ( |8.43 1 
in two places: first, explicitly in the square brackets; second, through {Dni\i)j_. Note, 
however, that the modulation of {Dnn)± does not produce any correction to the equilib- 
rium distribution function: ^-/^('^iv-^)/^ j^ ^ stationary solution of the equation no matter 
what happens to {Dnn)±- Thus, to find the linear in F^ correction to /IK/at), we should 
only include the modulation of con in the square brackets Eq. (8.43 1. As a result, 6fJl{lN) 
can be found from the equation 



d 

iw -I- — - {Dnn)j 
oIn 



_d_ 



+ l5{a)N-n) 



SMlN)=fiFa 



d {DNN)±e 



~li{a>N-li)lN 



dlff j3ia)N - /^) 



(8.58) 

Noting that the right-hand side of this equation is proportional to the eigenfunction mi(x), 
we easily find Sfjl and the susceptibility: 



Xn{(^) = 



1 



1 



1 - ico/r /3{ajN - ii)^ ' 



(8.59) 



Comparing this with the expression for S Ni(^), Eq. ( 8.54 1, we obtain the relation 



Sn{oj) ^ -z- \rciXN{(^), 
which is precisely the classical form of the fluctuation-dissipation theorem |[82]| . 



(8.60) 



8.4. Coupling to remote oscillators 

Here we analyze the perturbation theory terms which determine the diffusion in ac- 
tion of an oscillator which is sufficiently far away from the chaotic spot. Let the right- 
most oscillator of the guiding resonance be on the site n - 0, and let us focus on the 
oscillator at n = L, assuming it to be further away than the oscillators of the layer res- 
onance, i. e. L » Lf. As we have seen in Sec. |8.3[ relaxation of this oscillator is 
determined by the diagonal element Dn of the diffusion coefficient. The latter, in turn, 
is given by Eq. ( |8.28a i where the second term can be neglected as its denominator is 
determined by oscillators with \n\ ~ L( <^ L, and thus is large. 

To describe the coupling between the oscillators at « = and n = L for large enough 
L, it is sufficient to consider perturbation terms whose spatial structure is described by 

80 



diagrams of the kind shown in Fig. 12 First, consider Fig. 12 a). Namely, A^,- - 1 
nonlinear vertices are placed on some sites n'j, . . . , «^ ^ between and L, and each of 
these sites has m^*^, = 2(-l)^, while m^ = l,m''^ = (-1)^' , and the rest of m^ = 0. Indeed, 

k 

Vr <xi T^ which is clearly the lowest possible order in r. The number of resonances grows 

as 

(L-1)! n^le'^OiNJD^N, 



n= ^ kJ — '—[ , A^,-«L. (8.61a) 

(A^,-1)!(L-A^,)! \lnL^\ Nr j 

Next order in t is obtained by putting two oscillators at neighboring sites near one of 
the Nr- I nonlinear vertices, so it introduces a factor ~ A^^- in the number of resonances. 
Higher orders in r can be obtained if one separates the two ends of the solid lines near the 
^th nonlinear vertex by a distance l^. Assuming Ik «; L/Nr, which is equivalent to L^ - 
L <^ L,v/e can represent the number of resonances corresponding to such configurations 
as 

iNr-iy.(L-N,.y., Z. J"'+-+'~'-"^'-^ 

/i,...,//v,-i>o (8.61b) 

(L-1)! {Lr-L + Nr-2)\ 



{Nr-\y.{L-Nr)\{Lr-Ly.{Nr-2)\ 



This expression does not take into account changing the direction of the solid lines, 
shown in Fig. [T2]^b), which can be done at those nonlinear vertices where Ik > 0. Thus, 
Eq. ( 8.61b| ) gives a lower bound for the number of resonances, while the upper bound 



is obtained by multiplying it by 2^' for the two possible directions of each solid line. 
The number of terms of the perturbation theory which have the same order in r, p and 
contribute to the same resonance, is given by the product N = (h + l) . . . (In^-i + 1)- Note 
that Eqs. ( 8.61a i, (8.6 lb i do not fit in the discussion of Sec. 6.6 as an extra constraint is 
imposed (two sites of the resonance, « = and n = L, are fixed). 

It is not always necessary to invoke terms shown in Fig. [T2jb), as they bring an 



additional smallness r^(^'^^) to (Dll), while the resonances counted in Eq. (8.6 la I 
which correspond to Fig. [T2ja) may already be sufficient to remove the smallness of 
the Melnikov-Amold exponential, especially for very large L. Thus, at L > Lf one can 
write the following bounds: 

^ w e-^^ i-(i/^"P)r2(^-^^) (^? + . . . + ^^) < (D,,) < ^ r^^ (^f + . . . + ^l) . (8.62) 



Here each ^k is a product of random denominators, discussed in Appendix C The upper 
bound is obtained by simply setting L,- = L,Nr - 0, and the Melnikov-Arnold exponen- 
tial to unity. The lower bound is written by noting that if the optimal resonance is sought 
not on the full segment of the length L, but on a shorter one of the length L[, it will be 
similar to the layer resonance, so the corresponding factors are the same as those deter- 
mining the chaotic fraction, w, except that instead of the layer resonance one should take 
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(a) 

n=0 



n'2 



n=L 






n's 



(b) 

n=0 



n=L 





Figure 12: Examples of diagrams describing the coupling between the chaotic spot and a remote oscillator 
with A',. = 4. 



the next best drivi ng r esonance, hence the factor e~'^3'" (i/^'W ^j^j^ C3 ~ 1, as discussed 
in the end of Sec. 7.6 Also, w contains the thermal exponential e"^?^^, while Dn does 
not. Again, the disorder average of {Dn) diverges, as well as all higher moments, but as 
we will see in Sec. |9| it is the average of (Dil)'^ that enters the final result, so the sum 
^^ + . . . + ^?, can be effectively replaced by (const ~ 1)^. Moreover, even keeping this 
factor along with T^,p^ would be beyond our precision, so we set ^^ + . . . + ^?, ^ 1. 



8.5. Migration of chaotic spots 

The discussion of Sec. |8.2| regarded a segment of the chain containing one chaotic 
spot "living" on a given guiding resonance (m^, oj) = 0. The system dynamics in the 
space of the actions {/„) corresponded to the diffusion on the hypersurface determined 
by the constraints 



Z^" 



const. 



n 

(m^,a)) ^ 2j f^ii^n + gin) = 0. 



(8.63a) 
(8.63b) 
(8.63c) 



The last constraint is due to the fact that the system cannot leave the stochastic layer 



around the guiding resonance, and the variables p^, introduced in Sec. 8.2 are simply 
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(m«; to)=0 



Figure 13: A schematic representation of an intersection between two resonant hypersurfaces (irf, ai) = 0, 
(m*', w) = 0, shown by the smooth curves (the plane of the figure corresponds to the curved manifold 
/,o, = const, H,„, = const). The wiggly trajectories represent the system difi"using along the resonances. 
First, the system difi"uses along the g resonance and arrives at the intersection point. After some time it 
leaves the intersection point randomly choosing whether to continue along the g resonance or to switch to 
the g' resonance. 

the coordinates on the corresponding hyperplane. In the course of diffusion the system 
may cross another hypersurface, determined by the same action and energy conservation 
constraints Ifot = const, Htot = const, but with a different {irf ,a)) = 0. The time 
required for this to happen can be estimated as ~ (r/|/i|)^(Dii>~^e^s"'^. The prefactor 
corresponds to the time needed to cross the typical thermal region of the action space, 
and Dll is the diffusion coefficient of the most remote oscillator with mf,' t^ 0. This 



time coincides with the thermalization time of this oscillator, obtained in Sec. |8.3| The 
Arrhenius exponential corresponds to the number of attempts needed to overcome the 
activation barrier for the oscillators of the g' resonance. Thus estimate is valid for lengths 
L » Eg' IT; in the opposite limiting case the amount of action needed to put the system 
on the g' resonance is greater than the typical action, ~ LT l\ii\, stored in a segment of 
the length L, so the activation time is determined by the diffusion coefficient at lengths 

Once the system variables fall in the stochastic layer of the second resonance, we 
have two chaotic spots not far from each other, each one acting as a stochastic pumpp] 
Now nothing prevents either of them to leave its stochastic layer under the action of the 
other. Moreover, this happens in a relatively short time as the system only has to diffuse 



" In connection with the problem of system sticking to various structures inside the stochastic layer, 
discussed in the end of Sec. |5.2[ we note that when the system is changing the guiding resonance, sticking 
should be less probable than in the situation of the difi^usion along a single guiding resonance. Indeed, when 
two chaotic spots are present relatively close to each other, each of them effectively represents an external 
source of noise for the other one, preventing it from sticking. 
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across the thin stochastic layer. Once this happens, one of the chaotic spots is turned 
off, the other one remaining alone, and thus confined to its own stochastic layer. If it 
is the g' resonance that is turned off, there is no overall effect on the system. However, 
if it is the original g resonance, the system has effectively switched from one resonance 



to the other, as shown schematically in Fig. 13 As the two resonances are, generally 
speaking, located in different regions of the chain, this process corresponds to migration 
of the chaotic spot along the chain. 

In connection with migration of the chaotic spots it is worth mentioning a recent 
mathematical work 151]. For a one-dimensional chain of pendula with weak local cou- 
pling, energy above that of the unstable equilibrium was shown to be able to perform an 
arbitrary random walk along the chain. 

9. Current through breaks 

9.1. General remarks on the procedure 

Global transport properties of random one-dimensional systems tend to be deter- 
mined by rare strong obstacles ||53]| . The reason for this is that in one dimension such 
obstacles cannot be bypassed, in contrast to higher dimensions. Such rare obstacles were 
called "blockades" [64J, "weak hnks" [65 1, or "breaks" lUllllSl, and here the latter term 
will be used. In the present problem, a break is a region of the chain where the chaotic 
fraction w assumes anomalously small values on many consecutive sites, i. e., a region 
of the chain rarely visited by chaotic spots. Equilibration to the left and to the right of 
the break is assumed to be much faster than across the break, so the regions to the left 
and to the right can be assumed to be equilibrated with different values of the chemical 
potential and temperature, ni, Ti = 1//3l and fiR, Tr = 1/^r, respectively. The action 
and energy currents through the break, J-^ , 7^, are defined as the amounts of action and 
energy transferred from left to right per unit time. In the linear response approximation 
they can be written as 






Rl' 



Pr-ISl 



(9.1) 
It is convenient to borrow the terminology from the electric circuit theory, as discussed 



in Sec. 4.4 the 2-column on the left-hand side will be simply called the current, the 2- 



column on the right-hand side will be called the voltage drop across the break. The 2x2 
matrix R~y^ corresponds to the break conductance, and its inverse, R},, to the resistance. 
It is determined by the realization of disorder forming the break. 

Quite generally, validity of the linear response theory and finiteness of the coefficient 
7?^' are immediate consequences of the diffusive character of the system dynamics in the 
space of actions. In the present work the origins of the linear response theory are most 
easily seen in Sec. |8.3[ where a perturbation of the distribution function of an oscillator 
from its equilibrium value (say, putting it at a different temperature or chemical potential) 
is shown to correspond to an eigenfunction of the diffusion equation with a finite decay 
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rate. The analysis to be presented below corresponds to a more complicated situation, 
but the diffusion equation which was presented in Sec. [8] will lead to a definite expression 



iorRj/. 



In a stationary situation the current is constant along the chain, while the voltage 
drop across a segment containing many breaks is the sum of the voltage drops across 



each break. Thus, the macroscopic conductivity cr entering Eq. (3.5 1 can be expressed 
as 



-J 



Rt rib db 



(9.2) 



where Ub db is the spatial density (or probability per unit length) of breaks of a given 
type b. Eq. (9.2 1 is nothing but the usual addition rule for resistors in series, as dis- 
cussed in Sec. |4.4[ Analogously to the problem of electronic hopping conduction in 



one-dimensional samples ll53ll64ll65ll6Tll66ll . there are two strong competing factors: 
breaks with very large resistances are extremely rare, while frequent breaks have too 
small resistance. As a result, the dominant contribution to the integral comes from breaks 
which are close to a certain optimal configuration. This also means that the macroscopic 
resistivity cr~^ is self-averaging on the length scale corresponding to the typical distance 
between such optimal breaks. This section is dedicated to finding such optimal breaks 



and to evaluation of Eq. (9.2 1 analogously to Ref. |61 



In Eq. (9.2 1, Ri, is a random quantity, determined by the frequencies of all oscil- 



lators in the break region, while Ub db is the probability measure in the space of these 
frequencies. However, evaluation of Rb directly in terms of the frequencies seems to 
be hopelessly complicated. Instead, a break will be characterized by the corresponding 
spatial profile of the chaotic fraction |w„ = e~'^"}, defined by Eq. (7.8 1. Its probability 
distribution on a given site is given by p{A) = d[l - e~^^'^^]/dA, where S{A) was calcu- 
lated in Sec. |7] Assuming the chaotic fractions on different sites to be uncorrected, one 
can write the probability measure per unit length as 



Jnbdb^lj n d[l-e-'<'"^]. 



(9.3) 



|n|<L/2 



where S(i) is given by Eq. ( |7.51a| ), ^(i) - CiAP'pe'^'^'^^''^^^/^'f'^\ and L is a length which 
is much larger than the size of a break, but much smaller than the distance between 
breaks. The fact that the integral in right-hand side of Eq. ( |9.3| ) is proportional to the 
length follows from the translational invariance of the probability distribution and is 



discussed in detail in Appendix G 



Strictly speaking, Rb is not a functional of {A„} only, since A is determined by the 
guiding and layer resonances, while Rb depends on the diffusion coefficients, which are 
determined by the same guiding resonance as A, but also by the driving resonance which 
is different from the layer resonance. Still, it is quite strongly correlated with A, as 
discussed in Sec. |8.4[ Eq. (8.62i. In the following, the relation between Rb and {An} 
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Figure 14: The profile {A„] corresponding to some arbitrary break (shown by circles at integer positions). 
The solid curve shows the optimal break b^(x), corresponding to this profile [A,,], which touches [A,,] at two 
points (denoted here by «i and 112), and has the same resistance as the break with the profile IA„], up to 
power-law prefactors. 



will be approximated by a deterministic dependence Rh = /?({/}„}), and then the integral 



in Eq. (9.2 1 will be writen as an integral over {A„}. Determination of the dependence 
Rb = R{{An]) is the main task of Secs. |9.2| and |9.3| 

Obviously, the profile A„ corresponding to the optimal break must be symmetric with 
respect to its center, which can be assumed to be either at n = 0, or half-way between 
n = and « = 1 , without loss of generality. Also, at the center the profile A„ must have 
a maximum, /Iq » 1. Two mechanisms for the current flow through the break can be 
identified, which require different treatment. First, a chaotic spot may find itself near the 
center of the break, at \n\ ~ L™^ ~ Aq/ ln(l /r^'p) (the typical size of the layer resonance 
at the center of the break, as estimated in Sec. [7]l, which occurs very rarely, once in a 
time oc 6"^° . Second, a chaotic spot may reside somewhere on the slope of the break, 
\n\ ^ /}()/ ln(l /t''p) which occurs more frequently. However, the current will still be 
determined by the diffusion of oscillators near the center of the break, ~ Dqo, which has 
an additional smallness ~ t"^'"'. For the optimal break all these contributions must have 
the same order; indeed, if one of the contributions to the current is too small, this means 
that in the corresponding region of the break An can be decreased with no effect on the 
overall current, but with a gain in the probability. 

Putting the abovesaid in more formal terms, we introduce a one-parametric family 
of functions b;i{x), concave and symmetric, b^ix) = b;i{-x), which we will call optimal 
breaks, parametrized by their value in the maximum, ^i(O) = A. The resistance of a 
break with an arbitrary profile {A^} is the same (up to power-law prefactors) as that of a 
break from the optimal family, A"f' = b^{n), with the value of A determined as the largest 
A for which exists a shift xq such that /i„ ^ b;i{n - xq) for all n. Indeed, this choice of 
A implies that there are at least two sites, n = n\,n2, at which actually /i,,, = bx{n\), 

86 



bx{n2), and which lie on different sides of the break, as illustrated by Fig. 14 



Then the transport through such a break will be dominated by the chaotic spot when it 
resides on the sites n = n\,n2. The resistance of the break becomes a function of a single 



number A only, so that the integral in Eq. (9.2 ) can be written as 



<T 



' = f dAR{X)J{A)Y[e 



-s(b,m) 



(9.4) 



where the Jacobian 3{X) is analyzed in Appendix G The considerations of the previous 



paragraph already allow one to guess that the optimal breaks are likely to be triangular, 
Z?^(x) ~ A- 2\x\ ln(l/T), at least at large enough x, and thus the size of the break L}, ~ 
AI\n{\lT). 

9.2. Chaotic spot on the top of the break 



The definition of the currents given in Sec. 9. 1 leads to the following expression: 



j^H 



H>0 



1 



d{In)f 

dt 



+ 



H>0 





8l2 



dt 




0(r) 



(9.5) 



Here (...)/ denotes the average with respect to the distribution function /, whose time 



dependence is governed by the Arnold diffusion equation, Eq. ( 8.27 1. Integrating by 
parts and introducing the vectors iJ^ such that m^^q - 1, m^^q - 0, and a^ such that 



'fi _ 



oj ,. = oj„, w ,n - 0, we can write the current as 

n>\} "' n^U ' 



J' 



-2J 

kJk' 






Dkk' ^ W, W dpk. 



(9.6) 



where the diffusion coefficient D^t is given by Eq. ( 8.28a ), the distribution function / is 



defined in Eq. ( 8.26 1, and the stochastic layer width Ws is given by Eq. ( 7. 10a ). Note that 
it would not be correct to determine the distribution function itself from the stationary 
solution of Eq. ( |8.27 ). Indeed, the latter was written for a given fixed guiding resonance, 
while the position of the chaotic spot on the top of the break, considered here, is in fact 
very rarely occupied for short periods of time. The time required for the chaotic spot to 
jump off the top is of the same order as the equilibration time, so the distribution function 
on the left and on the right of the break is determined by a certain average over different 



guiding resonances [see Eq. (9.8 1 below]. 

Here we set / oc exp[-ySL//;(a»„ - //L/i?)^ - PLjRglnl^^ ^^^ '^^e oscillators to the 
left/right of the top of the break which do not participate in the guiding resonance and 
for which /?„ = /„, while the contribution of the oscillators participating in the guiding 
resonance will be simply neglected. The relative error introduced by neglecting these 
oscillators is small, since the number of the oscillators in the guiding resonance, 2Ng, 
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is much smaller than that for the layer resonance, 2Nt, and other resonances driving 
the diffusion (for these we expect A^,- ~ N[)- More important is the error introduced by 
setting the distribution function for other oscillators to its equilibrium value "by hands": 
it is likely to overestimate the current by a factor ~ 1, which is not so important as 



our main concern is about exponential factors. Using constraints (8.30 1, we express the 
currents as 



J' 



z 



1 



OJ. 



"1 



^nym LSR(a)„2 -/ir) -fidoJm - I^l)] 



(9.7) 



Here the thermal average (...) can already be performed over the equilibrium distribution 
where the difference between Pr,iir and PL,f^L is neglected, as we are working in the 
linear response regime. Moreover, the diffusion coefficient is always determined by 
the chaotic spot which has a specific position, so contributions from different possible 
positions n of the chaotic spots should be added with the corresponding weights, which 
we write as 



D 



n\n2 



= Z' 



,-^ n(") 



D 



« 1 (12 ■ 



(9.8) 



Li, the top 



As discussed in Sec. |9.1[ for the optimal break all terms of this sum should be of the 
same order. Since the dependence of D,",,^ on n is weak when \n\, \n\\, \n2\ 
of the break is smooth: bx{x) = Ao + 0(1) for \x\ ~ LJ"^" ~ Aol \n{llTPp). 
As mentioned in Sec. 



9.1 



^(«) 



the diffusion coefficient D„j„2 is strongly correlated with 
the chaotic fraction Wn- Indeed, both are determined by the same guiding resonance. 
The difference is that Wn contains the thermal exponential e~^s^^, while dJ,"^^ does not. 
Another difference is that while the chaotic fraction is determined by the best separatrix- 
destroying resonance, the diffusion coefficient is determined by the second best one, 
the criteria for the optimization being exactly the same. Indeed, the power spectrum of 



the stochastic pump, Eq. (8.6 1 contains the same Melnikov-Arnold exponentials as the 



stochastic layer width, Eq. (V.lOai. The difference between contributions of the best 



layer resonance and the second best resonance was discussed in Sec. 7.6 see Eq. (7.56 1. 
Thus, we adopt the following expression: 



D 



f.n) 
«1«2 



A3 



m„^m,^ 



-C2A,-C3\n^(llrPp) 



(9.9) 



with < C2 < 1, C3 ~ 1. Combining this with Eqs. ( |977) , (|9^, we obtain 

A3 



R 



-i| 



I top 



Z' 



-(l+C2M„-C3ln'(l/r"p) 






[ {nf,u^) {nf,o/) ]\. (9.10) 



It will be seen in Sec. 9.4 that for a typical break, contributing to the integral in Eq. (9.2 1, 
the chaotic fractions w„ are of the same order as for the optimal break, Wn ~ e"^'* ^"^ . This 



means that the sum in Eq. (9. 10 1 contains many (~ L™^") terms. Since m'''s are different 



for different n, the resulting matrix is not degenerate. Moreover, for each n the diagonal 
elements of the matrix can be estimated as {ifi,nf)^ ~ L™^", {{nf ,oJ^)) ~ L™^^/S^, and 
they are positive. The off-diagonal elements, {{a^ ,m[){nf ,i^)) ~ L^'"'A, have random 
signs. Thus, the summation over n effectively multiplies the diagonal elements by the 
number of the terms in the sum, ~ L™^", while the off-diagonal elements are multiplied 

by ~ /l™'"'. As a result, we arrive at the following expression for the contribution from 
top of the break: 



R- 



A' 



Kl+C2)^-C3ln^(l/r''p) 



I top 



1 



(L- 



max\-l/2 A 

A2 



(9.11) 



where the overall prefactor ~ (L™^")^ is omitted, as it is beyond our precision. In con- 



trast, the factor 
ments, is kept. 



(L™^") ^^^ determining the relative magnitude of different matrix cle- 



ft J. Chaotic spot on the slopes of the break 

When the chaotic spot resides on sites with \n\ > L™ ~ Ao/ ln(l /r^'p), its coupling 
to oscillators which are near the center of the break falls in the situation discussed in 
Sec. 8.4 The crucial point is that even for two neighboring oscillators the diffusion 
coefficient of the one which is closer to the chaotic spot is much larger than the diffusion 
coefficient of the more remote oscillator. This holds even if one takes into account that 
the diffusion coefficient is a strongly fluctuating random quantity whose moments are 
divergently 

Let us first consider a break which is symmetric with respect to the bond between the 
sites n = and n = I. Then the oscillator at « = is stronger coupled to the left part of 
the chain than to the right, while the oscillator at « = 1 is in the opposite situation. Then 
their distribution functions can be approximated by /o oc exp[-ySL(wo -fJ^Dh -/^Lgt^J^], 



f\ oc exp[-pR(a)Q- ij.r)Ii -/^Rglf /2], and calculation of the current according to Eq. (9.7| 
gives 
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I slope 
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(9.12) 



We have neglected the off-diagonal components Dio,Dqi of the diffusion coefficient, 
since the probability that an optimal driving resonance involves both these sites is low 
for \n\ » L^^ 



'-It is quite easy to check that for two independent random variables 1 < (1,(2 < °°, distributed according 
to a power law, p(() = a Q{( - l)/f '*", the probability of r^2 > (\ is given by t"/2. 



Let us now consider a break which is symmetric with respect to the site n = 0. Then 
the oscillator at n = is coupled equally well to the left and the right parts of the chain, 
so its distribution function /o(/o) should be found from the diffusion equation which has 



the same form as Eq. ( 8.43 1, but with two terms on the right-hand side, corresponding 
to the contributions from the left and from the right. In contrast to Sec. |9.2[ here this 
procedure is justified since the dynamics of the chaotic spots on the slopes of the break 
is faster than the dynamics of the oscillator in the center, so we can introduce the average 
diffusion coefficients 



D,^ 2] e-HD^^X Dr= Z ^'''-<<o>«- 



(9.13) 



,j<_^max 



Here the thermal averaging is performed with respect to the variables of oscillators with 
n < and « > 0, so that D^ and Dr depend on /q (namely, D^ r oc Iq). For the distribution 
function we obtain 



/o « exp 



DlPl + DrPr 



(^o/o + |/o') + 



Dl/3lI^l + DrPr/jr 



h 



(9.14) 



Dl + Dr \ "" 2 V Dl + Dr 

Denoting the average with respect to this distribution function by (. . .)o, and calculating 



the current according to Eq. (9.6 1, we write the conductance as 
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(9.15) 
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The two terms corresponding to the diffusion of the oscillators at n = ±1 are fully 



analogous to Eq. (9.12i. They are of higher order in t than the first term, but they 
are kept since the matrix appearing in the first term is almost degenerate. Indeed, for 
^UR K ^0 we can explicitly calculate the average: 



1 0)0 
0)0 0)0 

o»n = oio + 



/o 



IgT 



1 



0»n 



■,27^2 



w,', {i^'^Y + 2g^T'/icoo - liY + Oig') 



(/o>o , 



27.2 



(9.16) 



oJo-fi (o»o - nY 



+ 0{g% 



Note that the same kind of average actually appears in all terms in Eqs. (9.12 1, (9.15 1. 

If the size of the slopes is of the same order as the size of the top, ~ Lf^^, when the 
contribution from Eq. ( |9.12 ) or Eq. (9.15 1 is added to that from Eq. ( 9.11| ), they are of 
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the same order, so the resulting matrix is not degenerate, and can be inverted without 



problem to obtain R. In the opposite case, when the contribution from Eq. (9.12i or 



Eq. (9.15 1 dominates, matrix inversion should be done with care. 



Let us first consider Eq. (9. 12 1. The determinant of the resulting matrix is 



(DooXDiiXo;:, -ca\r + 2g^T'{(Doo) + (D„» 



(Dn) 



(Doo) 



(wo - fj-r (wi - fiy 



where we denoted 



<Doo>^ Y, e-'^{D^^), {Dn}= J] .-^"(Dlf), 



(9.17) 



n>L'^""' 



n<-L'p^ 



At first glance, the determinant can become small when ojq ~ a>i. Note, however, that 
both (Doo) and (Dn) contain factors l/(a»i - a>o)^, since coupling of the oscillator at 
« = 1 to the oscillators on the left slope necessarily involves tunneling through the 
oscillator at « = 0, and vice versa. Thus, already the first term in the determinant is 
regular, so the second one, which is small as 0{p^), can be safely neglected. 



For Eq. (9.15 1 the situation is somewhat different. The determinant is given by 



(Doo>(D„>(w'-w;)2 + (Doo>(D^ 



.i>K-^li)2 + (Doo>' 



(wo - ^ly 



where we defined (Dqo) = {DiDrI{Di + Dr))q, and (Dn), (D_i _i> are defined anal- 
ogously to Eq. (9.171. The first two terms have a smallness 0{t^) coming from {D\\), 
(D_i _i), while the last term has a smallness 0{p^). Still, the smallness of the determi- 



nant amounts to an overall power-law prefactor which is beyond our accuracy. 

9.4. Optimization 

For the moment, let us ignore the matrix structure of the break conductance R~^. 



Indeed, as we have seen in Sees. 9.2[ 9.3[ the matrix elements differ by a factor which is 
at most a power of T,p, or A, and possible problems with the matrix inversion amount to 
the same kind of factors. The optimization of the break shape involves exponentials, so 
it can be performed without paying attention to the matrix structure. 



From the discussion of Sees. 9.2 and 9.3 and the estimate (8.62i one can conclude 



that the optimal break functions b;i{x), defined in the end of Sec. |9.1[ should have the 



form, shown in Fig. 15 



bxix) = 



A + 0(1), 



\x\-AI\n{\lTPp), 



A - 2|x|ln(l/T), \x\ » AI\n{llTPp). 



(9.18) 



Then the product in Eq. \9.A\ is expressed in terms of [see Eq. (7.5 la I] 

bA{n) 



Y^S{b,{n)) = CipJ][bA(n)r exp 



Cln\l/TPp) 
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- C[pAP' exp 



Cln\l/TPp) 

(9.19) 
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Figure 15: The optimal break fc,i(x) at A » 1. bi(x) is defined as the family of functions, parametrized by 
A, such that the resistance of a break with an arbitrary profile [A,,] is the same (up to power-law prefactors) 
as that of an optimal break, A"/' = b^(n), with the value of A determined as the largest A for which exists a 
shift xo such that A„ ^ b^(n - xq) for all n (see Sec.|9.1|of the main text for details). 



where C\ is another logarithmic function of r and p. The main dependence of /? on /I 
is exponential, R{A) oc e'^i'^, 1 < Cj < 2 [see Eqs. (8.62i and (9.11 1], so the integral in 
Eq. (|9.4|) takes the form 



cr oc 



oo 

I dAexpic2 



A + lnJ{A)-C\pAP'exp 



Cln^il /rPp) 



(9.20) 



(9.21) 



The Jacobian is calculated in [Appendix G and it turns out that 

In J{A) ~ lnS{A) j^ , 

Cln\l/TPp) 

so it can be neglected in comparison with CjA. Then the optimal value /!*, maximizing 
the argument of the exponential, is found from the equation 



c-l^ + 



1 



A* ' C In^ {I /rPp) 
so that with logarithmic precision 



C[pAi ' exp 



Af 



Cln\l/TPp) 



7 1 1 

^ ^Cln^ In-. 

rPp p 



This determines the leading exponential in the macroscopic conductivity, 
Incr:^ =C'A 



C'^Cln^il/rPp) 



./LH. ^ w C'.Chi^ In -, 

l + (pi/A,)Cln^il/TPp) TPp p 



(9.22) 



(9.23) 



(9.24) 
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as well as the typical distance between such breaks, 

2 ;_; ^2. kC'C In^ — 

\ + {pilA,)C\rL^{\lTPp) ^ rPp 



InL, ? -^ ~ C'C In^ . (9.25) 



This distance diverges as t ^ or p ^ 0. Segments of the chain shorter than L* cannot 
be characterized by a conductivity, and the macroscopic equations given in Sec.[3]are not 
valid. For such segments only the total conductance of the segment can be defined as the 
response of the currents to the difference of potentials yu/T, l/T, applied at the ends of 
the segment. This conductance is determined by the strongest break on the segment, and 
thus is strongly dependent on the realization of disorder on the segment. Its probability 
distribution can be calculated analogously to Ref. {6V\, but is beyond the scope of the 
present work. 

The integration over A around the optimal value A^ can be performed in the Gaussian 
approximation, and the values of A effectively contributing to the integral are 

|^-^|~ln^. (9.26) 

It is easy to check that the cubic terms in the expansion of the argument of the exponential 
in /I - /Ih. are small as 1/ ln(l /r^p), so the Gaussian approximation is valid. Going back 
to the integrals over each individual An, we find that the integral converges at 

\An-b;^Xn)\-l. (9.27) 



This means that for a typical break, contributing to the integral in Eq. (9.2 1, the chaotic 
fractions w„ are of the same order as for the optimal break, w,, ~ g-^-uC") Hence, the 
conductance of such a break is contributed by many possible positions of the chaotic 
spot. 

Now one can turn to the study of the matrix structure of the macroscopic conduc- 
tivity (T. It also becomes clear when each of the two situation mentioned in the end of 



Sec. 9.3 is realized: the size of the top of the break, ~ /t* / ln(l /r'^p), and the size of the 
slopes, ~ /Ih./ ln(l/T), are of the same order unless ln(l/p) » ln(l/T), in which case the 
slopes are much longer. The inverse of a linear combination of matrices appearing in 



Sees. 9.2 and 9.3 always has positive quantities on the diagonal (~ A^ and ~ 1), while 
the off-diagonal elements contain terms proportional to a»„ plus a thermal correction, 
0{gT/{a>„ - jj)). This inverse should be averaged over the frequencies, as the macro- 
scopic conductivity is determined by many breaks. As a result, the diagonal elements 
remain of the same order of magnitude as before the averaging, while in the off-diagonal 
ones, ujn, which has a random sign, averages to zero, and only the thermal correction 
survives. Thus, the matrix structure of the conductivity is 



a oz 



1 0{gTIU) 

0{gT|\^A) 0(a2) 
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(9.28) 



We remind the reader that the matrix elements of cr have the following meaning [see 



Eq. (3.5 1]: cru determines the response of the action cun^ent to the gradient of fi/T 
when T is constant, crn determines the response of the action current to the gradient 
of 1/r when fi/T is constant, while cr2i and 0-22 correspond to the same for the energy 
current. 
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Appendix A. Thermodynamics of a strongly disordered nonlinear Schrodinger chain 

Consider a chain of a long but finite length L. Any microscopic state of the chain 
(a set of complex numbers {(A«)) can be characterized by definite values of the total en- 
ergy H, given by Eq. (23 1, and the total action (norm) Itot - H,, |iAnP (or their densities, 
'H = H/L and I = Itot/L, as introduced in Sec. |3]l. In a thermal state, characterized 
by two parameters, yS = l/T and fx, the microscopic variables are distributed according 
to the grand canonical distribution Pi{In,(f>n}) '^ g-/3(^-A'/w»)_ The temperature T and the 
chemical potential // determine the average action and energy per oscillator: 



I(T,fi) = - — In I exp -/3//({/„,0„))+^/i2]/„ ]~[ 
mT,fi)-fiI{M,T)-- — In I exp -/3H([I„,cf>n}) + /Sfx^^In ]~[ 

L n=l J n=l 



dl „ d(pn 
L 1 L 



(A. la) 
dl„ d(/)„ 



In 



(A. lb) 



Here the integration limits are ^ /„ < 00, ^ 0„ < In. In the limit of small r (strong 
disorder), Eqs. (A.la), ( A.lb| ) simplify: 



^^J Zn(r,yU) 

J 00 



^ ■l}(oj,-ii)I„-liglll2 ^dln ^ ^^^^^ 







Zn{T,fi) 



(A.2a) 



(A.2b) 
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r/(w„ - //), 



w„ - )U » gT 



(A.2c) 
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As usual, thermal fluctuations in the total energy and action are proportional to y[L, and 
are smaller than the average values themselves, LI and LfH, for a sufficiently large L. 
Moreover, even though each term in the sum is still a random quantity due to the ran- 
domness of ojn, summation over a large number of sites performs the eff'ective disorder 
average. Thus, in the above equations one can replace 



L ^1"^ 

1 \r-\ r d(jiJ„ 

"-^ -A/2 



(A.3) 



In other words, the total action and energy are self-averaging, and one does not have to 
invoke many disorder realizations. 

Performing the integrations, and keeping only the leading and subleading terms in 
the small parameter gT/{\fi\A) (i. e., expanding the exponential in g), we obtain 



r^^ A/2-// 

A "-A/2-// (/|2-A2/4)2 



'^-^^'^' .Oig^).Oir), 



(A.4a) 



These expressions work only for // < -A/2, |/i -i- A/2| » ^/gT (this is the condition 
for smallness of the subleading term with respect to the leading one). The assumption 
17^1 < (A/2)J corresponds to /j. not being too close to -A/2, and the strong inequality 



\1-(\ «: I A, Eq. (2.14i, translates into \fi\ » A. Expanding Eqs. (A.4a), (A.4b) in the 



small parameter A/\fj.\, we obtain 
I//I \ ^ 12//2 " 



^^^ ' + Oig^T^/fi^A^) + 0(AV/) + 0(T), (A.5a) 



fi^ 



'^ = ^ ^^^.o 2^ + O^g'l^^A^) + 0(AV/) + 0{T), 



(A.5b) 



which leads to Eq. (3.3 1. For // > -A/2 it is essential to keep glfj2 in the exponential. 



In this case instead of Eqs. ( A.4a i, ( A.4b ) one obtains 

r Ou + A/2)2/(2gA), -A/2 < // ^ A/2, 
1 f^/g, y" ^ A/2, 

, (2//3 + 3//2A/2 - AV8)2/(6gA), -A/2 < // =^ A/2, 
^ (ju^-A^n2)/2g, iu^A/2. 



(A. 6a) 
(A.6b) 



So far we were concerned the mapping (T,//) -^ (I, "H). It turns out that the inverse 
correspondence, (J, "K) — > (T,ii), is not always well-defined, as was noted in Ref. [311 
for the case without disorder. Below we study this problem for the strongly disordered 
case. 
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Figure A.16: The (J > 0,'7Y) half-plane in the dimensionless variables p = gl/A, h = g'H/A-. The solid 
line corresponds to T = 0; the region below this line is forbidden. The dashed parabola corresponds to 
r — > cx); the region above the parabola corresponds to non-thermal states which cannot be described by the 
grand canonical distribution. 



First, let us determine the allowed values of I and "K. For each J > the Hamilto- 
nian H has a minimum. Minimization of H under the constraint 2^ /„ = IL is performed 
by introducing the Lagrange multiplier //, and is equivalent to elimination of fi from 



Eqs. (A.6ai, (A.6bi: 



•H^ 



(J/6)( V32^ - 3A), ^ J ^ A/{2g), 
gl^ll - A\2Ag), I ^ Milg). 



(A.7) 



Clearly, there is no upper bound on "K for a given I. Indeed, one can consider a state 
h = IL, h = ... = II = 0, for which H = gI^L^/2 + 0{L), that is for L ^ oo the 
average energy density 'H can be made arbitrarily large. 

However, for thermal states the largest value of "K is obtained when T — > oo. In 
order for the partition function to remain finite one should keep yS/i = const < 0, then 
all In - I = -l/(ySyu), and we obtain the T ^ oo line to be TY = gl^ (in fact, it is the 
same as in Ref . ||58]| ; indeed, /3 = effectively removes the disorder from the partition 
function). Thus, the whole region 'H > gl^ of the right (J, "K) half-plane corresponds 
to non-thermal states (i. e., no T and fi corresponding to given 1,'H can be found in 



this region); see the discussion in Ref. ||58]| for details and Fig. A.16 for the graphical 



representation of the right (J, "K) half -plane. The assumptions Eqs. ( 2. 14 1, (2. 15 I ensure 
the fact that the considerations of the present work are restricted to thermal states only. 
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Appendix B. On the spreading of a finite-norm wave packet 



Here we are going to see what the macroscopic equations of Sec.[3]predict when an 
initial condition in the form of a finite-size wave packet is taken for I{x), 'Hix). Let us 
use the variables p, h, introduced in Eq. (|3.2|), and assume 



(B.l) 



h <0, p «c \h\ «: p. 
Then we can rewrite the transport equations in a closed compact form: 



d 



d{tA) 



P 
h 



d_ 
dx 



-C'ln^(l/r''p)ln(l/p) 



1 

Cip 



Ci 



d_ 
dx 



-\Ip 

hliUp^) 



Cl,C2~l. 



(B.2) 

Now one can estimate the growth of the typical size of the cloud, L{t), by writing p ~ 
poLo/L{t) and approximating the spatial derivatives as dx -^ i^/L{t) -^ p/(poLo). Let us 
also assume ln(l/T) <sc ln(l/p). Then one obtains the following estimate: 



dp 



-1 



dt 



-C'ln^p- 



L{t) 



,[(C'r' In ?]"3 



(B.3) 



which is Eq. (3.8 1. Here we rescaled the time as t 
integral 



th.l{p^L^) and approximated the 



/ 



e^'i">Jp 



-1 



.C'lnV- 



as power-law prefactors are beyond our precision. 

It is important to check that if inequalities ( |B.1| ) are satisfied at the initial moment, 
they remain satisfied for all times. This is certainly the case if the p cloud and the h cloud 
have the same length scale L{t) during the expansion: p ~ poLo/L{t), h ~ hoLo/L{t). Let 
us check that in the course of the expansion the two length scales Lp, L/, do not separate. 
Clearly, L/, » Lp is impossible since the expansion of the h cloud outside the p cloud 
is blocked by the exponential. To see that L/j «; Lp is also impossible, let us write the 
currents as (the terms proportional to Ci are small and can be neglected) 



JP 

J' 



(C2/U)d,h 



dxP 

- iC2/6)(h/p)dxP 



(B.4) 



Thus, the p cloud expands on its own, so one can study the expansion focusing on p only. 
If at some instant L/j «: Lp, then dh/dt will be larger than dp/dt by a factor LilL^, so 
the expansion of the h cloud will be faster, and it will eventually catch up. 



Note that the length scale L*, Eq. (3.7 1, determining the validity of the macroscopic 
transport equations, diverges much faster than L. This means that, strictly speaking, 
the spreading of a wave packet cannot be described by macroscopic equations. Still, 
as discussed in Sec. 4.5. 5[ the effect of the fluctuations at short length scales is just to 
modify the factor C' , so the functional form of the dependence L{t), Eq. (B.3 1, should 
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not be affected. However, the spreading of the wave packet is subject to strong sample- 
dependent fluctuations. 

For not sufficiently large L another point is worth checking. One of the key assump- 
tions of the present work is that the system can overcome an arbitrarily high activation 
barrier E, in a sufficient time, oc e^^^ , i. e., an arbitrarily large amount of action and 
energy can be concentrated at a given point at some instant of time. However, the wave 
packet contains only a fixed finite amount of action and energy, so sufficiently large fluc- 
tuations are forbidden. The condition for applicability of the Gibbs distribution can be 

written as 

H-iil.„. . 1 1 

(B.5) 



H-JUltot „, . T 1 



» C In^ ^- In 

I ' rPp p 

The right-hand side of this inequality represents the typical value of the activation barrier 



corresponding to the optimal break, as discussed in Sec. 9.4 The left-hand side is just 
the size of the cloud, L. Since p oc 1 /L, this inequality is satisfied after a sufficiently 
long time of the expansion. 



Appendix C. Probability distribution of a product of denominators 

Here we study the properties of the probability distribution 



^niO = J 



SU 



1 



Xi 



dx] 



. ClAjfi — 



^2 („-!)!■ 



(C.l) 



If one adopts Stirling's formula for (n - 1)!, Eq. (C.l I can be rewritten as 



PniO 



^y/2nin- 1) 



exp 



2{n - 1) "^ 3{n - 1)2 



(C.2) 



When all terms in the exponent except the first one are much smaller than unity, which 
is the case if 



In 



^ 



at-l 



« 3^'\n - 1) 



2/3 



(C.3) 



they can be neglected, and Eq. (C.2i corresponds to the log-normal distribution. The 
total probability contained in the log-normal region is determined by the error function 
erf([(9/8)(« - 1)]'^^). Thus, the rest of the probability decays as e~" , i. e., slower 
than exponentially. The average of ^^, which diverges for the actual distribution, for the 
log-normal one is finite and equal to e^^"~^\ 

Another interesting region is the tail of the distribution, where the probability integral 
can be calculated by the steepest descent: 



CO OG 



-l+{n-l)\nl_ 



dl 
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In"^ 



iB^ 



(«-!)! ^(«-l)!ln(^/e"-i) 
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ln^(^/e«-i) 



-I- 



(C.4) 



where we have linearized the argument of the exponential with respect to I. The correc- 
tion is small when 



In 



^ 



,n-\ 



» 



V«-i. 



(C.5) 



Comparing this with Eq. (C.3 1, we see that the two regions actually overlap. The approx- 



imate form of the distribution PniO ^^ the tail region can be obtained by differentiating 



Eq. (C.4 1 with respect to ^. 

Let us now study the sum of N independent random variables ^ with random signs, 
z = I ± ^1 ± . . . ± ^yvl- As all moments of the distribution PniO diverge, the central 
limit theorem is not valid and the distribution of z, p^iz), does not tend to a Gaussian. 
In the tail region (large z), Pziz) coincides with the distribution of max{^i, . . . ,^/\f}, so 
the sum is simply determined by its largest term. The maximum of p^{z) is at z = 
0, due to the possibility of the cancellation between different terms in the sum, even 
though each individual term is bounded from below [p„{0 vanishes for ^ < 1]. These 
properties will be discussed in more detail in the end of this appendix. Here we note 
that the random variable max{^i, . . . ,^^} represents an upper bound on z since their 
distributions coincide at large z, while the probability of small values of max{^i , . . . , ^yvl 
is suppressed, in contrast to small values of z, as discussed above. 

The probability for max{^i, . . . , ^yv) to be smaller than some value z is given by the 
product of the probabilities. 



'P{^u.-.,(n<z} 



-\pn 



\N 



iOd^ 



exp 



N 



ln"z 



z(«-l)!ln(z/e«-i) 



(C.6) 



Most of the probability is contained in the region of z where the argument of the expo- 
nential is of the order of unity, i. e. 



1 



ln"z 



1 



z(n-l)!ln(z/e"-i) A^ 
At large n the solution of this equation is given by 

-1-1/ In /(//'/") 



where the function f{x) is defined as the solution of the equation 

/ 



1+ln/ 



= X. 



(C.7) 



(C.8) 



(C.9) 



In particular, if N ~ c", then z ~ {ef{c))". To give an idea of the behavior of /(x), we 
note that /(I) = 1, and x\n(ex) ^ f(x) < 2A;ln(ex) f or x > 1. In the region around 



the solution of Eq. (C.7i the probability ( |C.6[ ) can be approximated by linearizing the 
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logarithm of the argument of the exponential in Eq. (C.6l with respect to Inz, as in 



Eq. (C.4i: 



'P{^u...,^N<e- 



exp 



-exp 



In/ 
1+ln/ 



{A + n\n{ef)) 



(CIO) 



A lower bound on z can be obtained if one replaces Pziz) by a box distribution 
Pz{Q)6{\IPz{Q) - z) thus suppressing the probability of large z- Unable to calculate 
^^(O) explicitly, we obtain the lower bound in three steps: (i) replace PniO by the cor- 
responding log-normal distribution, which already suppresses the probability of large z; 
(ii) apply the central limit theorem to the sum of log-normally distributed random vari- 
ables; (iii) replace the resulting Gaussian by a box of the same height. This procedure 
results in a box distribution with a height larger than /j^CO). To account for random signs 
of ^k, we redistribute the log-normal symmetrically between the positive and negative (, 
then the second moment of this distribution is given by the average of ^^ over the log- 
normal, e^^""'\ so the central limit theorem immediately gives 



Pz(0) < 



^^e2(n-l) 



(C.ll) 



rff-l 



P^iO 



1^1 > ^cx 



(C.12) 



To conclude this appendix, let us return to the distribution of the sum of random 
variables not satisfying the central limit theorem. Consider a random variable ^ whose 
probability distribution p^iO is symmetric, p^{0 = P({-0^ ^nd for |^| > ^oo can be 
approximated by 

where ^o ^ ^oo so that the total probability contained in the region \x\ > ^oo does not 
exceed unity. Let 1 < a ^ 3 so that the total probability converges, but the second 
moment diverges. Then the central limit theorem is not applicable for the sum z = 
^1 -I- ... -I- ^yv- The probability distribution for this sum can be expressed as 

oo r oo "I-'V oo r oo 

Pz(z) - J^^''' Je-''^pa)d^ « J^e'^'exp -N J{\ - e'''^) p^{Od^ ■ 

— OO 1:^00 J —OO L —OO 

(C.13) 
At small ^ «; 1 /^co the ^ integral can be approximated as 



- cos k^ 



Je-''^p,iOd^ = 1 - C' / ^-^^^ d^ + OikYj = 



^a 



1 - (^^o) 



a- 1 



oo 

J- 



- COSM 



du + OCfc^-'^C) 
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+ Oik'^i) 



(C.14) 



The integral over u is equal to -r(l - a) sm{nal2). The 0{k'^^^) term can be neglected 
when N » (^oo/^o)*^^~"^''*^""^\ so we can write 



Pziz) 



CO 

/ 



dk 



Izl » Zn, 



2n~ i (W2Ko"Vz", 

where zn represents the typical value of |z|: 



(C.15) 



Zn - (o 



-Nr{\ - a) sin 



na 



l/(a-l) 



(C.16) 



The expression for p^{z) at large \z\ » Zn was obtained by expanding the exponential 
as e"''^^'^'" = 1 - \kzN\"~^ + . . ., performing analytical continuation from the positive 
imaginary semiaxis of z, and using the relation T{a) r( 1 - a) sin na = n. It turns out to 
coincide with the tail of the probability distribution of max{|^i|, . . . , l^yvll (up to a factor 
of 2, which accounts for the difference between the distributions of z and |z|): 



tt-r 



!„,<,,.. |e.p-^|. 



N^r' 



M-\ 



(C.17) 



Let us now consider the sum of positive random variables, assuming 7?^ (^) to vanish at 
^ < and to have the form ( C.12| ) without the factor 1/2 at large ^. For 1 < a < 2 
the integration is analogous to Eq. (C.14i, but the result must be analytical in the lower 



complex half -plane of k. Then, instead of Eq. (C. 15 1 we obtain 

QO CO 

dk 



Pz(z) 



! 



Jkz+Nr{l-a)iik^of 



In 



- r — e-"^ Imexp [-;Vr(l - a) (/<^o)''"^e'™] , 



(C.18) 
which at large z again gives NC^'"^ Iz" ■ For 2 < a < 3 the first moment, ^, of /7f(^), is 
finite, so it should be subtracted, and at large z we obtain NQ~^ l{z - Nl)". 

Appendix D. Change of variables in the vicinity of a resonance 

Consider N oscillators, n = I,. . .,N, and suppose that the condition 

N 



fl=l 



dH(h,...,iN) 

din 



(m, o)) - 0. 



(D.l) 



is satisfied for some values of the actions / = /''. Then the slow phase is introduced as 



2_^^n'Pn- 



(D.2a) 
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To complete the transformation to new phases 0[, we write 

N 

01 - 0> 'P'k>l = /^Uknlpn, (D.2b) 

n=\ 

where the matrix u^n is chosen so that the vectors m^ = {u^i, . . . , uun) satisfy the orthog- 
onality relations 

i^k>vU^'>\) = ^kk', iiik>vm) = 0. (D.3) 

Then, defining u^ = m, we can easily invert the matrix: 

I Ukn^ K > i, ^—f 

^ n=\ 

and write the completeness relation as 

N 



ninm, 
\m 



-^rr + > , UknUkn' - 5nn'- (D-5) 



J:=2 



The transformation for phases fixes the transformation to the new actions pk in the 
vicinity of the resonance point (so that the whole transformation is canonical): 



N 

L = YjP'^^k- (D-6) 



k=i 

This gives 

dH ^^ dH din ^ OH 



"Pk ^ din dpk ^ din 

so that dH/dpi = 0, and pi is the effective pendulum momentum. The kinetic part of the 
pendulum Hamiltonian is then determined by d^Hldp\ - g\m]^, and this is what appears 
in Eq. ([73]). 

Let us now approximate 

q2t7\"l / a2i 



dH \^\ { d ti 
dpk [dpj) \dpkdpi 



^^(^) = ^^ + lirTl I^^Tl^PUO, (D.8) 



where the derivatives are taken at the point /''. The second term appears when the equa- 
tions of motion are written to the first order in pi which gives d^[/dt = {d^Hldp\)p\. 
Inverting Eqs. (D.2ai, (D.2b i, we obtain 

^ , dH I d^H Y^ d{m,oj) . 

<t>nit) = Y,Ut)ukn = ^J^[m,^^,m]^ -^m, (D.9) 
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which is Eq. ([0]l. 

Let us now see how the stochastic layer thickness translates into an effective integra- 
tion measure on the resonant manifold. Namely, we consider the integral of a function 
F{Ii, . . ., /„), smooth on the scale of the stochastic layer width: 

J n -^^ J ^7^ ^ " J ^^^('"I'^i + ■ ■ • + mNcbN)Y\dI,„ (D.IO) 
*=^ layer "=1 

where At is the measure to be determined. Since the transformation is canonical, 

k=l n=i 

M is determined from the condition 

r dpi d(p' r 

W.S = — = dp]_M6(miaJi + . . . +mNCJN), (D.12) 

layer 

where Ws is the volume of the stochastic layer in the pendulum phase space. In order to 
resolve the 5-function we write 

d{m,cjj)^ } , m„———dln' ^\m,—— m\dpi+ } \u,,—— m\dpk (D.13) 

The last term vanishes, and the integration over p\ gives 

M = W, y mn 4-^ mn' = g\m\^W,. (D.14) 



It is this measure that appears in Eq. ( |7.8| ). 

Appendix E. Integrals used in thermal averaging 

Let us consider the following integral: 

CO 

Jai...aM) = I exp -xi - . . . - x„ - -^^ — \dx[... dx,,. (E. 1) 



At A » 1 the leading exponential asymptotics is easily found by noticing that the maxi- 
mum of the exponent in the integrand is reached at 

/ . \l/(ai+...+Qf„+l) 

4'"' = ^A-, A , (E.2) 

\a ' ...an") 
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so we obtain 

In /«!...«„ (A » 1) - -(ai + . . . + «„ + 1) 
Let us note that 



a] a 



l/(ai+...+Q-„ + l) 



(E.3) 



(ai + . . . + a„) In '— ^ \ o-yt In a/t ^ (ai + . . . + &„) ln(ai + . . . + a„). (E.4) 

k=i 



n 



The left inequahty is the Jensen's inequality following from the convexity of the function 
xlnx (the equality holds when all ak = a), while the right inequality obviously follows 
from the monotony of In x (the equality holds when all a^ except one are zero). Denoting 
ai + . . . + a„ = na, we obtain an estimate 



exp 



a- 



na + \ 

na/(na+l) 



a/{na+l) 



^ Jai...aM » 1) ^ exp 



na + \ 



(j^Q,\na/(na+l) 



1 l/(na+l) 



(E.5) 



Even though the saddle point approximation is not supposed to be valid when some 
of the ats are zero, the right inequality is still correct since the integration over the 
corresponding n - \ coordinates simply gives j e'^^'^dxt = 1, while the integration over 
the remaining coordinate in the saddle point approximation gives the same result, as can 
be checked directly. 

Let us now focus on the case when when all a^ = a, 



Ja{A) ^ I exp 








Ji, j^ ... -^n 



(Xi . . . XnT 



Ct^l ... CI,J\.iJ^ 



(E.6) 



and improve this estimate by calculating the prefactors in front of the exponentials in the 
Gaussian approximation: Let us perform an orthogonal change of variables 



n-l 



(E.7) 



where the unit basis vectors e'" are orthogonal between themselves, as well as to the 
direction of y\\ : 



Z' 



k'^k 



■ Z' 







(E.8) 



k=\ k=i 

To the second order in the perpendicular deviations. 



Xl . . . Xfi 



l/3'll 



3^11 



^Jnj 2 \ ^Jn 



n~2 f!-l f 



\k+k' 
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,m'=\ ^ 
-1 

2^3'L 



mm 
^k^k' 



m,m'=l 
n-2 n-\ 



Jmjn 



(E.9) 



where we took advantage of the orthogonaUty: 

( m \ ( m \ 



z 



mm 



U=l 



\k'=\ 



E^^ Z^^' -Y^^yk--^. 



k=\ 



Then the integral is easily calculated as (F = y\\l ^/n) 



Ja(A » 1) - f ^fiidY fd''-^y^exp 



mm' 



,-f |2 



-nY- 



a A 13/xl 



2 yjiQ' y2 

^^)"^ " r e""^"'^''^"" y("«+2)("-i)/2 jy _ 
aA)"-i j 



?i(q'?i + 2 - a) / (27r)" 



(■Q,^yi/(na+l) g^p 



na + \ 



a 



{aA) 



\l(na+l) 



2 V na + 1 

(E.IO) 

The third line is obtained from the second one by denoting u - y(««+2)(«-i)/2+i^ ^^^ 
integrating over u in the Gaussian approximation. The prefactor is large, so to obtain a 
lower bound we can simply omit it. Then, setting n = 2Ng - m^s and a - 1/4, we arrive 
at Eq. ( |7T31 ). 

To obtain the upper bound, we note that ;c"^^e~'"''^ ^ e""''^, which gives 

Ja{A » 1) ^ n{na + 2){2nlef''^ exp 



{aA) 



1/(00'+ 1) 



(E.lla) 



Note that the integral over yj_ can be calculated in the Gaussian approximation only 
when aAj Y"" » 1 , otherwise it is of the order of unity. Thus, to obtain an upper bound, 
we just set the factor {Y"" I aA)"^^ in the second line of Eq. (E.IO I to unity and integrate 
over F" = m in the Gaussian approximation: 



/ (2nYrfi 
Ja{A » 1) ^ ^/ 1—1— -(aA)(2«- !)/(««+!) exp 



na -\- \ 
^ n{%nle^)"l'^ exp 



na + \ 



a 



■{aA) 



\l{na+l) 



< 



(E.llb) 



{aA) 



l/(na+l) 



where the last inequality was written using {xe '^^^)" ^ (2/e)". Since Stt/^^ - 3.4 . . . > 



27r/e = 2.3 . . ., we adopt the latter estimate. This is precisely Eq. (7.20 1 



Appendix F. Effective diffusion coefficient atong a thin tayer 

Consider the diffusion equation for the distribution function f{W, p2, ■ ■ ■ , Pn) 
f{W, p)msi flat space (the volume element being dWdp2 ■ ■ ■ dpN)'- 



df 
dt 



= -d'^J'w - djcJk, Jw - -Dwwdwf - Dwkdkf, Jk = -Dkwdwf - Dkk'dwf. 
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(F.la) 



Here d^ = d/dW, dk = d/dpk, and the summation over repeating indices ^ = 2, . . . , A^ is 
assumed. The diffusion is confined to a thin layer W\{p) < W < Wzip), W2ip) -Wi{p) = 
Wsip), so that the normal component of the current vanishes at the boundaries: 



JW Z ■'*: 

OPk 



0. 



(Rib) 



W=Wi2 



All dependencies on p/c are assumed to be slow on the scale Ws'. W^Okf «: /■ Thus, one 
can seek the distribution function in the form 



W2(p) 

Aw, p)^f{p) + Rw, p), J f{w, jT)dw = o. 



(F.2) 



Then our goal is to eliminate / and find the effective equation for f(p). 

Integrating Eq. (F. la) over W and using the boundary conditions, we obtain 

W2{p) 

W,(/?)^^ - -dk f JkiW,if)dW. 



(F.3) 



Let us perform the gradient expansion, / - /^'^ +/*^^^ + 0{dl). To deal with the diffusion 



|2\. 



equation and the boundary condition we need Jk to 0{dk) and Jy/ to 0{d^): 

-Jk = bkk'dk'f + DkwdwP^ + 0{dl), (F.4a) 

-Jw = Dwkdkf + DwwdwF^ + DwkdkF^ + DwwdwF^ + 0{dl). (F.4b) 



Let us multiply Eq. (F.lal by Ws and subtract from it Eq. (F.3 1. Since dt = 0{d^), we 
can neglect dtf and write 

= Wsdw [Dwkdkf + Dwkdkf + Dwwdwf] + 

( W2 



+ WsdkDkk'dk'f -dk 



! 

W2 



Dkt dW 



dk'f+ 



(F.5) 



+ WsdkDkwdwf 



-dkf 

Wi 



DkwdwfdW + Oidh. 



There are only two terms of the order 0{dk), from which we obtain 



/ 



^(1) _ 



( ^^ ^ 
Dwk 



fS^dW 
J Dww 



dkf, -Jk - I Dkk' - ^J^^ I dk-f - Dl^,dk-f. (F.6) 



D 



WW 
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Substituting Jt^ into Eq. (F.3 1, we obtain the desired effective equation for f{p): 



dt W,{p) 



W2(P) 



- r 

dpk J 



dW 



Dkk'iW,jf) 



wm 



Dkw{W,iT)Dwk'{W,{f) 



dfip) 
dpk' 



(F.7) 



Finally, we have to make sure that the correction /^^^ is regular (its explicit form is not 



even needed). Collecting the 0{d^) terms in Eq. (F.5 1, we obtain 



dwDwwdwf =dwDwkdk 



+ — dk 



The boundary conditions translate into 



Dwwdwf \w=w, 2 ~ ^^'"'^'^ 






dW 



dk'f-dkDl,,dk'f+ 



(F8) 



Wi 



( ^1.2 ^ 



dk'f. 



f^dW 



dW 



1,2 



dk'f+^r^Dl,dtf. 
opk 



(F.9) 



Eq. (F.8l is compatible with these boundary conditions and can be straightforwardly 
integrated, indeed giving f^'^^ of the second order of smallness. 

Appendix G. Jacobian for integration over breaks 

Our goal is to find the Jacobian J'{A) in the transformation 

-S(hAn)) 



fY[p^^n)dAn^ ^ dAJ{A)Y[ 

^ n ^ n 

where the probabihty distribution /j(/l„) is represented as 

d[l - £'-'^(''«^i 



p{A„) = 



dAn 



S\An)e 



-S(A„) 



(G.l) 



(G.2) 



The calculation is based on three main observations. 

First, as discussed in Sec. |9.1[ for each break profile {An} one can specify the two 
sites, «i,«2> where it touches the corresponding optimal break, /l„,j = b^ini^ - -^o). 
for some values of A,xq, which are determined by the profile {An}. Thus, the whole 
integration manifold of {An} can be spUt into submanifolds corresponding to given values 
ofn\,n2'. 

~ ~ "2) (G.3) 

Ml ^02 ' 



Jji^ J]~[^^^ 2] JjI("''«2 
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Second, at fixed ni,n2 there is one-to-one correspondence between /l,n,/l„, and the pa- 
rameters A, xq of the optimal break b^in - xq). That is, for an arbitrary profile {A,,} the 
parameters A, xq are in fact determined by four numbers n\,n2. Am , An^- The transforma- 
tion (/l„, , /l,,^) <r^ {A, xo) is determined by 

dArn^ = dAbximj - xq) dA - b\{ni2 - xq) dxQ, (G.4) 

where we denoted b'j^{x) = dbx{x)ldx. The integration measure transforms as 

dAmdAn^_ - \b\{ni - xq) 5^^i(«2 - -^o) - b\{n2 - xq) dAbxini - xq)] dxo dA. (G.5) 

Third, at fixed n\,n2,Arn,An2, or, equivalently, at fixed ni,n2,A,XQ, the integration range 
for An with n t n\,n2is simply bx{n-XQ) < /l„ < oo, so the integration is straightforward. 
Thus, one can write 



I cf^l^K^)^ /li I dAnidAn.^ p{Ani) p{Am) Y\ exp[-S(Zji(?i-xo))] = 

n 11 ^«2 ni^rn,n2 

2 f dAn,dAn,S\An,)S\An,)e'-^'\ 



Here we introduced 

^-s(i) ^ j-j ^^p [_5(^^(^ _ ^q))] ^ (G.7) 

assuming it to be independent of xq. This is not strictly true, since n is integer and xq 
is continuous, so it is actually a periodic function of xq with unit period, but we neglect 
this modulation. 

Let us now pass from A,^,A„2 to A,xo keeping A fixed, and change the summation 
variables as ni = n, n2 = n + I as well as the integration variable as xq = « -i- / - x: 

"2 

^ I dxoS '{bx{ni - xq)) S'{b^{n2 - xq)) X 

[b\{ni - Xo) dAbA{n2 - xq) - b\{n2 - xq) dAbxini - xo)] = 



" 1 

X 



oo /,„ 



^m 



(G.8) 



dx 



dS {b^ix - 0) dS (bAix)) dS {bAx)) dS {b,{x - 0) 



" ^ ^ J "■" dx dA dx dA 

n=-oo 1=0 J^ L 

- LJ(A), 

which gives precisely the sought Jacobian. Due to translational invariance, the sum over 
n produces the total length L, which appears in Eq. (9.3 1. The value of l,„ax - ImaxW is 
determined from the equation bx{lmaxl'^) - 0, and the range of x integration is given by 

, ? (G.9) 

, _ i-max ^ ^ . l-max ^ i ^ i 

' ~ 5^ ''C 5% •■maxt n ^ ' ^ I-max- 
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Replacing the summation over I by integration, we can rewrite the Jacobian as 

OS (bA-y)) dS (b^x)) dS {bAx)) OS (bA-y)) 







J(A) = - dxdy 



+ 



dy dA dx dA 



I 12 (G.IO) 



.„) I ^~.., 



J 

" ''max I ^ 



References 



[1 

[2; 

[3 

[4: 

[5 
[6 

[7 

[8: 

[9 
[10 

[11 
[12: 

[13 

[14: 

[15 
[16 

[17: 

[18 
[19 
[20 
[21 

[22: 
[23: 

[24 
[25 

[26: 
[2?: 
[28: 
[29: 

[30 
[31 
[32 
[33 
[34 
[35 
[36: 
[37: 



P. W. Anderson, Absence of diffusion in certain random lattices, Pfiys. Rev. 109 (1958) 1492. 

A. Lagendijk, B. van Tiggelen, D. A. Wiersma, Phys. Today 62 (2009) 24. 

M. E. Gertsenshtein, V. B. Vasil'ev, Theory of Probability and its Applications 4 (1959) 391. 

N. F. Mott, W. D. Twose, Adv. Phys. 10 (1961) 107. 

J. Frohlich, T Spencer, C. E. Wayne, J. Stat. Phys. 42 (1986) 247. 

S. A. Gredeskul, Y. S. Kivshar, Phys. Rep. 216 (1992) 1. 

D. L. Shepelyansky, Phys. Rev. Lett. 70 (1993) 1787. 

M. 1. Molina, G. P Tsironis, Phys. Rev. Lett. 73 (1994) 464-467. 

M. I. Molina, Phys. Rev. B 58 (1998) 12547-12550. 

T. Paul, P Leboeuf, N. Pavloff, K. Richter, P Schlagheck, Phys. Rev. A 72 (2005) 063621. 

T. Paul, P Schlagheck, P Leboeuf, N. Pavloff, Phys. Rev. Lett 98 (2007) 210602. 

A. lomin, S. Fishman, Phys. Rev. E 76 (2007) 056607. 

J. Bourgain, W.-M. Wang, J. Eur. Math. Soc. 10 (2008) 1-45. 

S. Fishman, A. lomin, K. Mallick, Phys. Rev. E 78 (2008) 066605. 

G. Kopidakis, S. Komineas, S. Flach, S. Aubry, Phys. Rev. Lett. 100 (2008) 084103. 

A. S. Pikovsky, D. L. Shepelyansky, Phys. Rev. Lett. 100 (2008) 094101. 

S. Fishman, Y. Krivolapov, A. Soffer, J. Stat. Phys. 131 (2008) 843. 

M. Albert, T Paul, P Leboeuf, N. Pavloff, Phys. Rev. Lett. 100 (2008) 250405. 

S. Tietsche, A. Pikovsky, Europhys. Lett. 84 (2008) 10006. 

S. Flach, D. O. Krimer, C. Skokos, Phys. Rev. Lett. 102 (2009) 024101. 

C. Skokos, D. O. Krimer, S. Komineas, S. Flach, Phys. Rev. E 79 (2009) 056211. 

S. Fishman, Y Krivolapov, A. Soffer, Nonlinearity 22 (2009) 2861. 

A. lomin, Phys. Rev. E 80 (2010) 037601. 

H. Veksler, Y Krivolapov, S. Fishman, Phys. Rev. E 80 (2009) 037201. 

T Paul, M. Albert, P Schlagheck, P Leboeuf, N. Pavloff, Phys. Rev. A 80 (2009) 033615. 

M. Mulansky, K. Ahnert, A. Pikovsky, D. L. Shepelyansky, Phys. Rev. E 80 (2009) 056212. 

M. Albert, T. Paul, N. Pavloff, P Leboeuf, Phys. Rev. A 82 (2010) 01 1602(R). 

J. D. Bodyfelt, T. Kottos, B. Shapiro, Phys. Rev. Lett. 104 (2010) 164102. 

H. Veksler, Y Krivolapov, S. Fishman, Phys. Rev. E 81 (2010) 017201. 

Y. Krivolapov, S. Fishman, A. Soffer, New J. Phys. 12 (2010) 063035. 

A. lomin, Phys. Rev. E 81 (2010) 017601. 

C. Skokos, S. Flach, Phys. Rev. E 82 (2010) 016208. 

S. Flach, Chem. Phys. 375 (2010) 548. 

M. Mulansky, A. Pikovsky, Europhys. Lett. 90 (2010) 10015. 

T V. Laptyeva, J. D. Bodyfelt, D. O. Krimer, C. Skokos, S. Flach, Europhys. Lett. 91 (2010) 30001. 

A. Pikovsky, S. Fishman, Phys. Rev. E 83 (2011) 025201. 

M. Johansson, G. Kopidakis, S. Aubry, Europhys. Lett. 91 (2010) 50001. 

109 



[38] T. Schwartz, G. Bartal, S. Fishman, M. Segev, Nature 446 (2007) 52. 

[39] Y. Lahini et al, Phys. Rev. Lett. 100 (2008) 013906. 

[40] Y. Lahini et al, Phys. Rev. Lett. 103 (2009) 013901. 

[41] J. Billy et al, Nature 453 (2008) 891. 

[42] G. Road et al. Nature 453 (2008) 895. 

[43] B. Deissler et al. Nature Physics 6 (2010) 354. 

[44] A. Dhar, J. L. Lebowitz, Phys. Rev. Lett. 100 (2008) 134301. 

[45] A. Dhar, K. Saito, Phys. Rev. E 78 (2008) 061 136. 

[46] M. Mulansky, K. Ahnert, A. Pikovsky, Phys. Rev. E 83 (201 1) 026205. 

[47] V. Oganesyan, A. Pal, D. A. Huse, Phys. Rev. B 80 (2009) 115104. 

[48] I. Garcia-Mata, D. L. Shepelyansky, Eur. Phys. J. B 71 (2009) 121. 

[49] D. O. Krimer, R. Khomeriki, S. Flach, Phys. Rev. E 80 (2009) 036201. 

[50] W.-M. Wang, Z. Zhang, J. Stat. Phys. 134 (2009) 953-968. 

[51] V. Kaloshin, M. Levi, M. Saprykina, (to be published). 

[52] N. R Mott, Philos. Mag. 19 (1969) 835. 

[53] J. Kurkijarvi, Phys. Rev. B 8 (1973) 922. 

[54] B. V. Chirikov, Phys. Rep. 52 (1979) 263. 

[55] A. J. Lichtenberg, M. A. Lieberman, Regular and Stochastic Motion, Springer- Verlag, New York, 

1983. 

[56] D. J. Thouless, J. Phys. C 5 (1972) 77. 

[57] N. Bilas, N. Pavloff, Eur. Phys. J. D 40 (2006) 387. 

[58] K. O. Rasmussen, T. Cretegny, P. G. Kevredikis, N. Gr0nbech-Jensen, Phys. Rev. Lett. 84 (2000) 

3740. 

[59] G. M. Zaslavsky, Chaos in dynamic systems, Harwood Academic Publishers, Newark, 1985. 

[60] B. V. Chirikov, V. V. Vecheslavov, J. Stat. Phys. 71 (1993) 243. 

[61] M. E. Raikh, I. M. Ruzin, Sov. Phys. JETP 68 (1989) 642. 

[62] V. K. Melnikov, Soviet Mathematics-Doklady 3 (1962) 109. 

[63] V. I. Arnold, Soviet Mathematics-Doklady 5 (1964) 581. 

[64] V. K. S. Shante, C. M. Varma, A. N. Bloch, Phys. Rev. B 8 (1973) 4885. 

[65] P A. Lee, Phys. Rev. Lett. 53 (1984) 2042. 

[66] A. S. Rodin, M. M. Fogler, Phys. Rev. B 80 (2009) 155435. 

[67] M. Khodas, S. Fishman, O. Agam, Phys. Rev. E 62 (2000) 4769. 

[68] A. N. Kolmogorov, Dokl. Akad. Nauk SSSR 98 (1954) 525. 

[69] J. Moser, Nachr. Akad. Wiss. Gottingen, Math.-Phys. Kl. 1 (1962) 1. 

[70] V. I. Arnold, Uspehy Math. Nauk 18 (1963) 13. 

[71] C. Froeschle, J.-P Scheidecker, Phys. Rev. A 12 (1975) 2137. 

[72] E. Diana, L. Galgani, M. Casartelli, G. Casati, A. Scotti, Theor. Math. Phys. 29 (1976) 1022. 

[73] M. Casartelh, E. Diana, L. Galgani, A. Scotti, Phys. Rev. A 13 (1976) 1921. 

[74] J. Poschel, Commun. Math. Phys. 127 (1990) 351. 

[75] D. M. Leitner, P G. Wolynes, Phys. Rev. Lett. 79 (1997) 55. 

[76] V. Y. Demikhovskii, F M. Izrailev, A. I. Malyshev, Phys. Rev. Lett. 88 (2002) 154101. 

[77] I. L. Aleiner, D. L. Altshuler, G. V. Shlyapnikov, Nature Physics 6 (2010) 900. 

[78] J. D. Bodyfelt, T. V. Laptyeva, C. Skokos, D. O. Krimer, S. Flach, Phys. Rev. E 84 (2011) 016205. 

[79] G. Czycholl, B. Kramer, A. MacKinnon, Z. Phys. B 43 (1981) 5. 

[80] M. Kappus, R Wegner, Z. Phys. B 45 (1981) 15. 

[81] A. B. Rechester, R. B. White, Phys. Rev. Lett. 45 (1980) 851. 

[82] L. D. Landau, E. M. Lifshitz, Statistical Physics, Part 1, Pergamon Press, Oxford, 1959. 

[83] H. B. Callen, T. A. Welton, Phys. Rev. 81 (1951) 34. 



110 



